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EXECUTIVE  SUMMARY 


This  report  describes  the  results  of  a  study  of  several 
numerical  methods  for  calculating  points  on  the  distribution  of  a  sum 
of  statistically  independent  random  variables.  The  report  is  directed 
to  practitioners  of  statistical  and  numerical  methods.  Inunediate 
motivation  for  the  study  arose  in  connection  with  the  random  time  to 
accomplish  a  collection  of  tasks.  However,  application  to  a  variety 
of  problems  is  anticipated  because  of  the  generality  of  the  metnods. 

Study  objectives  are  to  investigate  the  relative  accuracy  and 
computational  effort,  viz,  run  time,  for  eacn  of  the  following  methods: 

(a)  Evaluation  of  closed-form  solutions  for  particular  cases. 

(b)  Discrete  numerical  convolution  of  probability  densities. 

(c)  Normal  probability  approximation  to  the  distribution. 

(d)  Numerical  Inversion  of  the  Laplace  transform  of  the  convolu¬ 
tion.  (Bellman's  method). 

(e)  Erlang  approximation  for  convolutions  of  a  two-parameter 
WeibuU  distribution.  (Johnson's  method). 

(f)  Convolution  of  probability  densities  using  the  FFT  algorithm 
for  calculating  finite  Fourier  transforms. 

(g)  Monte-Carlo  simulation. 

Normal  approximation  for  sams  of  Independent  random  variables 
(RVs)  is  made  in  several  areas,  including  quality  control  and 
analytic  network  theory.  Because  of  the  frequently  uncritical 
assanption  of  Normality,  the  error  of  this  approximation  Is  a 
particular  focus  here. 

Methods  are  sketched  for  deriving  analytic  expressions  for 
the  distribution  of  the  sura  of  RVs  of  certain  distributions.  Each 
numerical  method  is  described  and  Illustrated  using  RVs  from  several 
distributional  forms,  such  as  uniform,  exponential,  gamma,  and 
Wei  bull,  as  well  as  mixture  models.  In  terras  of  run  time  and 
accuracy,  some  methods  are  particularly  suited  to  certain  distribu¬ 
tional  forms.  If  problem  applications  are  quite  special  and  If 
the  time  for  program  ceding  (as  well  as  running)  is  a  consideration, 
Monte-Carlo  simulation  may  be  tne  preferred  method.  All  computer 
source  programs  are  listed  In  annexes. 
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SUBJECT:  Methods  for  Calculating  the  Probabil  v  Distribution  of 
Suras  of  Independent  Random  Variables 
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2.  Background 


Problems  associated  with  the  distribution  of  sums  of  independent 
random  variables  (RVs)  occur  in  various  analyses.  Examples  are: 

(a)  estimating  the  total  cost  of  an  end  item  or  a  project,  given 
component  cost  estimates;  (b)  estimating  time  to  complete  a  series 
of  sequential  tasks,  (This  problem  can  be  generalized  to  estimating 
the  completion  time  for  a  series  of  networks.);  (c)  estimating 
dimensional  variability  in  an  assembly  of  serially  arranged  parts; 
and  (d)  estimating  statistical  confidence  limits  on  the  mean  of  a 
random  variable  (RV).  Changing  problem  context  may  obscure  the  math¬ 
ematical  identity  of  these  familiar  problems.  Often,  the  probabil¬ 
ity  distribution  of  the  sum  is  assumed  to  be  Normal,  since  the 
central  limit  theorem  guarantees  Normality  as  the  nviaber  of  RVs 
in  the  sum  becomes  infinitely  large.  However,  if  either  tail  of 
tne  distribution  of  the  sum  of  a  small  number  of  independent  RVs 
Is  to  be  estimated  with  accuracy,  it  is  prudent  to  be  cautious  in 
Immediately  assuaing  Normality.  This  report  addresses  the  issue 
of  accuracy  of  a  Normal  approximation  and  other  issues  associated 
with  different  methods  of  calculating  the  camulative  distribution 
function  (c.d.f.)  of  a  sum  of  n  random  variaoles  (RVs),  when 
the  RVs  have  a  variety  of  distributional  forms. 


$ 


I 


3.  Study  Objectives 


Specific  objectives  of  the  study  reported  here  are:  (a)  identify 
the  error  of  approximation  for  the  c.d.f.  of  an  n-component  sura  is 
n  Increases;  This  error  is  examined  for  cases  In  which  all  of  the 
components  have  the  same  distribution  and  for  cases  in  which  tne 
form  of  the  distribution  is  the  same  but  in  which  tne  par;imeterj  are 
unique,  (b)  obtain  closed-form  expressions  for  the  c.d.f.  of  the 
sum  for  special  cases,  wnlch  may  ue  used  to  check  various  numerical 
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metnods;  (c)  obtain  measures  of  computational  effort  for  several 
numerical  methods  for  comparative  purposes;  and  (d)  suggest  which 
methods  are  suitable  for  treating  particular  cases. 

4.  Discussion 


The  sum  of  two  RVs  has  a  distribution  which  is  the  mathematical 
convolution  of  the  component  distributions.  Thus,  the  sum  of  n 
RVs  has  a  distribution  which  is  the  n-fold  [1]  convolution  of  the 
component  distributions.  The  problem  of  calculating  the  dist¬ 
ribution  of  the  sum  is,  then,  equivalent  to  obtaining  an  n-fold 
convolution.  Tnis  fact  can  be  exploited  in  calculating  the  c.d.f. 
of  the  sum  numerically  and  in  obtaining  a  formula  for  the  c.d.f.  of 
the  sum.  Consider  the  case  of  two  continuous  RVs,  defined  on  0  to 
Infinity.  Call  the  variables  x  and  y  and  call  the  sum  z.  The  c.d.f. 
of  any  RV  will  be  denoted  by  F(.),  with  a  subscript  referring  to  the 
variable  of  interest.  Similarly,  the  notation  for  the  probability 
density  function  (p.d.f.)  of  interest  will  be  f{.),  with  a  specifying 
subscript.  Thus,  tne  c.d.f.  and  p.d.f.  of  the  RV  x  are,  respectively, 

F  (x)  rnd  f  (x).  For  this  case,  the  convolution  theorem  yields: 

X  X 

F  (z)  5  lntegral(0,z):  F  (z  -  t)  f  (t)  dt  .  (1) 

z  y  X 

Similarly,  from  (1),  the  p.d.f.  of  z  is  written  as 

f  (z)  z  lntegral(0,z):  f  (z  -  t)  f  (t)  dt.  (2) 

z  y  X 


For  specific  distributional  forms  the  indicated  integration  may 
be  simple  to  carry  out.  If  so,  a  closed  form  expression  for  the 
desired  convolution  is  obtained.  If  not,  one  can  use  discrete 
numerical  convolution.  The  numerical  equation  is  obtained  from 
equation  (1)  by  discretizing  the  domains  of  the  functions  at,  say, 
m  identical  points:  t(i},  for  1  ie  i  le  m.  Then,  the  diferential 
form  f  (t)  dt  is  replaced  by  a  probability  difference  and  the 

X 

integral  becomes  a  sum,  as  follows: 

F  (z(k))  z  Sum  lzl,k:  F  (z(k)-t(l))  (F  (t(l)-F  (td-l))).  (3) 
z  y  X  X 

Tne  accuracy  of  the  numerical  method  improves  by  increasing  the 
numoer  (m)  of  discrete  points,  assuming  that  the  range  t(m)  -  t(l) 
adequately  covers  the  domain  of  tne  c.d.f.  of  z  in  the  sense  that 
tne  upper-tall  probability  beyond  z(m}  is  negligible— say,  1/100,000. 

If  numerical  convolutions  are  to  be  performed  recursively,  it  is 
necessary  to  anticipate  the  domain  of  the  hlghest-order  convolution 
wnen  choosing  z(m).  Because  of  the  need  for  a  high  density  of  discrete 
points,  me  size  of  a  generally  becomes  quite  large  for  four  or  more 
convolutions.  Clearly,  this  situation  produces  a  computational 
ourden  which  increases  rapidly  with  the  order  (n)  of  convolution. 

For  n  greater  tnan,  e.g.  4,  other  numerical  methods  may  be  preferred 
on  tne  basis  of  efficiency. 

"  order *of "convolution  is  defined  nert  as  the  number  of  distribu¬ 
tions  being  convolved.  This  is  equal  to  the  i  of  RVs  in  a  sum. 
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5.  Integral  Transforms 


In  dealing  with  convolutions  It  Is  helpful  to  use  a  theorem  [1]  from 
the  theory  of  Integral  transforms— either  Laplace  or  Fourier.  Tnat 
theorem  states  that  the  transform  of  a  convolution  of  two  function  is 
the  product  (in  the  complex  plane)  of  the  transforms  of  the  functions. 
Following  a  UK  convention,  I  denote  the  Laplace  transform  of  a  function 
with  the  function  symbol  having  an  asterisk  superscript.  For  example, 
the  Laplace  transform  of  f(x)  is  f*(s),  with  complex  argument  s. 

Thus,  the  p.d.f.  of  z  in  (2)  can  be  characterized  by  the  transform: 

f*(s)  =  f»(s)  f»(s)  .  (4) 

2  y  X 

If  the  RV  2  is  added  to  another  RV  w  yielding  the  sum  v,  one  can 
immediately  write  the  transform  of  the  p.d.f.  of  v  as 

f*(3)  =  f»(s)  f»(3)  f»(3)  ,  (b) 

V  y  X  w 

instead  of  convolving  f  (y)  with  f  (x),  and  the  result  with  f  (w). 

y  X  w 

Of  course,  it  is  necessary  to  be  able  to  invert  the  transform  to 
achieve  the  desired  result.  More  will  be  said  about  this  later. 

For  many  probabl'ity  distributions  of  interest,  the  Laplace  transform 
can  be  written  simple  form.  Examples  are  the  uniform  distribution 
on  (0,a),  which  has  the  Laplace  transform  of  the  p.d.f.: 

f*(3)  s  (1  -  exp(-as ) )/a/s,  (6) 

and  the  exponential  distribution  with  rate  parameter  r,  whose  p.d.f. 
transform  Is 

f«(s)  =  r/(s  ♦  r).  (7) 

If  each  of  the  n  random  variables  in  the  sum  has  the  same  distribution, 
the  transform  of  the  p.d.f.  of  the  sum  is  Just  the  ntn  power  of  tne 
transformed  p.d.f.  For  the  sum  of  n  uniform  (0,1)  deviates,  yielding 
tne  RV  t. 


n  n 

f»(s)  s  (1  -  exp(-3))  /s  .  (d) 

t 

The  Laplace  transform  of  tne  c.u.f.  of  t  is  obtained  from  the 
transformed  p.u.f.  simply  by  dividing  by  s,  since  tne  c.d.f.  is 
just  the  integral  of  the  p.d.f.  To  facilitate  inversion,  the 
expression  for  the  n  th  power  of  1  -  exp(-3)  Is  bo  expanded  as 
a  sum  of  blnjalal  terms: 

I 

Sum  over  i  {0,n);  C(n,l)  (-1)  exp(-is), 

where  C(n,i)  is  tne  #  of  c^xablnat ions  of  n  objects  taken  1  at  a  time. 

ffT  An Vxi^sl t ion  of  the  theorem  is  found,  e.g,,  in  denkins,  G.M.  and 
Watts,  D.d.  Spectral  Analysis,  Holden-Day,  c.  H03. 
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The  transformed  c.d.f.  of  t  can  be  readily  inverted  analytically: 

1  n 

F  (t)  =  Sum  over  i  (0,n):  (-1)  C(n,i)  u(t-i)  (t-i)  /n!  ,  (9) 

t 

where  u(t-x)  is  the  unit  step  function  at  x.  This  expression  is 
quickly  and  accurately  evaluated,  even  for  large  values  of  n  (n  >  10). 
Calculation  was  performed  by  the  routine  NFOLD.U  on  the  Prime  9955 
minicomputer  with  a  run  time  limited  by  the  print  buffer,  l.e.,  in  a 
fraction  of  a  second.  This  program  is  listed  in  Annex  A. 

6.  If  the  RV  of  interest  (t)  is  the  sum  of  n  identical  exponential 
RVs,  the  Laplace  transform  of  the  p.d.f.  of  t  is,  from  (7)$ 

n  n 

f«(s)  =  r  /(s  ♦  r)  .  (10) 

t 

This  expression  also  has  a  simple  inverse: 
n  n- 1 

f  (t)  s  r  t  exp(-rt)/(n-1)?  (11) 

t 

This  is  recognized  as  a  gamma  p.d.f.  with  shape  parameter  n  and 
rate  parameter  r.  This  result  illustrates  the  familiar  theorem 
that  the  sum  of  n  Identical  exponential  RVs  has  an  Erlang  distri¬ 
bution,  l.e.,  a  gamma  distribution  with  integer  shape  parameter. 

It  follows  immediately  from  (10)  that  the  sum  of  N  identical  gamma 
distributions,  having  shape  parameter  n,  is  also  a  gamma  distribu¬ 
tion  with  snapo  parameter  Nn,  since  the  Laplace  transform  of  its  p.d.f. 

Nn  Hn 

f*(s)  *  r  /(s  ♦  r)  ,  (12) 

t 

nas  the  same  form  as  the  transform  of  the  gamma  p.d.f.  in  (10). 

Aitho  the  Laplace  transforms  in  these  examples  have  simple  inverses, 
transforms  are  still  useful  in  calculating  convolutions  of  probability 
distributions  when  this  condition  does  not  exist.  The  reasons  for  this 
assertion  are:  (a)  that  the  transforms  of  the  distributions  being 
convolved  ire  often  simple  functions  of  s,  (b)  that  the  product  of 
such  transforms  are  easy  to  evaluate,  and  (c)  that  numerical  methods 
exist  for  calculating  the  inverse  Laplace  transform.  One  such  method 
was  developed  by  Ricnard  Bellman  (Ref  (la)).  I  have  found  this  method 
useful  in  several  applications,  such  as  In  solving  integral  equations 
(Ref  (1b])  as  well  as  for  obtaining  the  distribution  of  sums  of  RVs. 
Farmer  discussion  of  Bellman's  method  is  deferred  to  a  later  point. 

7.  Measures  of  Accuracy 


Several  measures  can  be  used  In  describing  the  accuracy  of  numerical 
methods  for  approximating  the  c.d.f.  of  a  sum  of  RVs.  Two  are  used 
here:  (a)  the  mixlmum  absolute  error  over  a  finite  set  on  the  domain 
of  the  c.d.f.,  and  (o)  the  square  root  of  the  mean  squared  error  or 
HiS  error,  evaluated  over  the  same  set  of  points.  For  most  methods 
In  this  study,  I  nave  used  20,  equally- spaced  points  on  the  domain  of 
th<*  c.d.f.,  such  that  tail  probabilities  are  *ess  than  0.01  beyond 
tnv  range  of  points  used.  An  exception  to  this  selection  of  points 
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is  made  when  using  Bellman's  method.  In  that  case  16,  log-spaced 
points  are  used  to  span  the  range  of  the  sum.  As  an  illustration 
of  these  measures,  consider  the  Normal  approximation  to  the  N-fold 
convolution  of  a  standard  (0,1)  uniform  distribution.  Using  the 
exact  result,  given  in  (9),  the  measures  of  error  are  calculated  for 
several  choices  of  N.  These  are  shown  in  Tabie  1. 

TABLE  1 

ERROR  IN  THE  NORMAL  APPROXIMATION  OF  CONVOLUTIONS  0?  N 
STANDARD  UNIFORM  PROBABILITY  DISTRIBUTIONS 


N  Max  Abs  Error  RMS  Error 


2'  ^oroffit  ■  ^0.00^ 

3  0.0097  0.0054 

4  0.0074  0.0038 

5  0.0057  0.0029 

10  0.0028  0.0013 


In  many  applications,  the  error  associated  with  a  very  large  (small) 
value  of  the  c.d.f.  is  more  appropriate  than  either  of  the  above 
error  measures.  For  5  convolutions  of  a  uniform  distribution,  the 
error  of  a  Normal  approximation  is  about  0.1<  for  values  of  the 
c.d.f.  >  0.95.  By  nearly  any  measure,  4  or  5  convolutions  of  a  given 
uniform  distribution  is  well  approximated  by  a  Normal  distribution 
whose  mean  and  variance  are  N  times  the  uniform  mean  and  variance. 
However,  not  all  distributions  of  sums  of  uniform  RVs  are  this  well 
approximated  by  a  Normal  c.d.f.  Tne  case  of  sums  of  different 
uniform  RVs  is  considered  below.  One  may  ask  if  Monte-Carlo  simul¬ 
ation  is  competitive  in  terms  of  accuracy— if  not  in  terms  of  run 
time — with  a  Normal  approximation.  For  the  case  considered  above, 

20  thousand  Monte-Carlo  replications  produces  a  typical  RMS  error 
of  0.002  to  0.003.  This  is  about  the  same  accuracy  as  the  Normal 
approximation  for  N  =  5.  The  run  time  for  20,000  replications  on 
the  Prime  9955  is  approximately  a  linear  function  of  the  number  of 
RVs  in  the  sum.  For  this  case,  approximate  run  time  T,  in  seconds, 
is  given  by 

T  =  8(N  -  2)  ♦  30.  (13) 

For  simple  cases  such  as  this,  Monte-Carlo  Is  quite  expensive  in 
terms  of  run  time.  However,  Monte  Carlo  becomes  more  attractive 
when  the  problem  becomes  mathematically  intractable. 


8.  Sums  of  Non- identical  Uniform  RVs 


The  n-foia  convolution  of  the  standard  uniform  distribution  was 
obtained  in  closed  form  (9)  by  inversion  of  tne  Laplace  transform, 

given  in  (9).  Tnls  result  ca^  ^  generalized  by  permitting  each 

of  the  n  anlform  distributions  to  have  a  different  range,  but  with 
common  threshold  parameter.  Thus,  the  ktn  member  of  the  set  is 
defined  on,  say,  0  to  3(k).  The  Laplace  transform  of  tne  p.d.f.  of 
the  sum  (t)  is,  then. 


f*(3)  =  Product  over  k=l  to  n:  (1  -  exp(-a(k)s )/ (a(k)3 ) .  (14) 

t 

5 


Tne  inverse  transform  is  somewhat  complicated  to  obtain,  and  is 
not  derived  here.  The  exact  c.d.f.  for  the  sum  of  n  different 
uniform  RVs  is  simply  presented,  with  the  following  definitions,  as 

n  k 

r  (t)  =  l/a/n![t  +  Sum  over  k  (1,n)!  (-1)  Sum  over  j  (1,C(n,k)) 

t 

n 

u(t  -  S  (n))(t  .  S  (n))  ]  ,  (15a) 

kj  kj 

where 

a  =  Product  over  k=l  to  n:  a(k) 
and  where 

S  (n)  is  the  j  th  sum  of  the  k  tuple  of  n  values  of  a(i),  taken 
kj 

k  at  a  time.  For  example, 

S  (n)  r  a(J),  1  le  j  le  n, 

S  (n)  r  a(1)+a(2)  and  S  =  a(n-l)+a(n).  (15b) 

21  2C(n,2) 

The  implementing  computer  program,  given  in  Annex  B,  calculates  the 
error  of  a  Normal  approximation  for  a  larger  class  of  sums  of  uni¬ 
form  RVs.  Consider  the  following  special  case  in  which  the  range 
of  the  k  th  uniform  RV  in  the  sum  of  n  is  k.  Normal  errors  in  the 
c.d.f.  of  the  sum  are  given  in  Table  2.  Compare  with  Table  1. 

Note  that  the  errors  are  about  twice  those  in  Table  1.  However, 
even  these  errors  are  relatively  small  (1X  or  less)  for  N  =  5. 

TABLE  2 

ERROR  IN  THE  NORMAL  APPROXIMATION  OF  CONVOLUTIONS  OF  N 
DIFFERENT  UNIFORM  PROBABILITY  DISTRIBUTIONS  » 


“n  Max  Abs  Error  RMsTrror 

’2 . •oTftTf - 'Offfff - 

3  0.0179  O.OlOil 

a  0.0131  0.0071 

5  0.0102  0.0054 

10  0.0049  0.0024 


•  Range  of  the  k  th  uniform  RV  is  taken  to  be  k. 
h.  Sums  of  Identical  Exponential  RVs 


For  a  somewnat  different  picture,  consider  the  case  of  a  sum  of  N 
exponential  RVs  from  the  same  c.d.f.  The  error  of  a  Normal  approxima 
tlon  Is  shown  in  Table  3  as  a  function  of  N. 


TABLE  3 

ERROR  IN  THE  NORMAL  APPROXIMATION  OF  N  CONVOLUTIONS  OF 
AN  EXPONENTIAL  PROBABILITY  DISTRIBUTION 


N 

Max  Abs  Error 

RMS  Error 

“2” 

TTotilff"' 

3 

0.0769 

0.0365 

4 

0.0698 

0.0316 

5 

0.0596 

0.0278 

10 

0.0916 

0.0183 

15 

0.0390 

0.0192 

Tne  errors  shown  in  Table  3  are  about  one  order  of  magnitude 
greater  than  those  in  Table  1,  indicating  that  sums  of  exponential 
RVs  approach  Normality  much  more  gradually  than  suras  of  uniform 
RVs.  For  a  sura  of  15  exponentials,  the  Normal  error  is  about  M 
or  less  for  values  of  the  c.d.f.  >  0.98.  Clearly,  this  example 
indicates  a  need  for  caution  in  applying  the  Normal  assumption. 

10.  Suras  of  Different  Exponentials 


The  p.d.f.  of  the  sura  of  exponential  RVs  from  the  same  distribution 

was  shown  (11)  to  have  the  Erlang  form.  If  a  set  of  n  exponential 

RVs  from  distributions  with  unique  mean  values  are  summed,  the  form 
of  the  c.d.f.  is  somewhat  complicated.  However,  an  analytic  model 
exists  for  this,  more  general  case.  The  computer  program  which  is 
used  for  evaluating  this  distribution  is  found  in  Annex  C.  If  the 
rate  parameter,  r(k),  of  the  distribution  of  an  arbitrary  kth  RV  Is 
unique,  the  c.d.f.  of  the  sura  of  n  RVs  is  given  by 

F  (t)  =  r  Sura  over  i  (1,n):  A(i)(l  -  exp(-r(l )t ) )/r(i ) ,  (16) 

t 

where  r  =  Product  over  k=1  to  n;  r(k), 

and  where  the  vector  A(*)  is  the  solution  of  a  certain  matrix 

equation:  M  A  =  B.  Elements  of  the  B  vector  are  all  zero  except 
the  nth  (last).  A  typical  element  of  M,  m(i,J),  involves  the  sura 
of  all  (i-1)  tuple  products  of  r(k),  with  k  not  s  to  J.  Thus,  e.g., 

ra(3»j)  =  Sum  over  k  ne  J  (l,n):  Sum  ever  1  >  k,  ne  j;  r(k)r(l). 

The  Ist  row  of  M  has  elements  =  1.  Other  rows  are  like  tne  one  above. 
Equation  (16)  can  be  used  to  calculate  the  Normal  c.d.f.  error  for  a 
special  case.  In  a  set  of  n  exponential  RVs,  let  the  range  of  the 
mean  values  be  fixed  at  2.  Let  the  kth  RV  have  the  mean  value 
1  (k-l)/(n-1).  Tne  Normal  errors  for  the  c.d.f.  of  tne  sura  of  these 

RVs  are  shown  in  Table  4.  Comparison  with  the  results  of  Table  3 
indicates  that  greater  errors  of  Normal  approximation  occur  whan 
the  RVs  in  the  sura  have  different  mean  values.  For  this  example 
the  error  is  about  10$  greater  than  for  the  n-fold  convolution  of 
the  same  exponential  distribution. 


TABLE  4 

ERROR  IN  THE  NORMAL  APPROXIMATION  OF  CONVOLUTIONS  OF 
N  DIFFERENT  EXPONENTIAL  PROBABILITY  DISTRIBUTIONS* 


Max  Abs  Error 


RMS  Error 


"O.IOBr 

0.0823 

0.0693 

0.0634 

0.0436 


0.0477 

0.0393 

0.0338 

0.0297 

0.0195 


kth’RV'fn  a  set  or~n,“et *the“niean  value “Ge'Oc-TTTTn- 1 ) . 


U.  An  Exponential  Mixture  Model 


In  forming  a  sum  of  independent  RVs,  one  may  think  of  each  RV  as 
the  duration  of  a  particular  activity  in  a  serial  network  of  n 
activities.  Tne  distribution  of  the  sum  is,  then,  the  distribu¬ 
tion  of  completion  time  for  the  network.  A  variation  of  this  model 
is  one  in  which  the  nodes,  separating  activities,  permit  two  exit 
patns,  each  of  which  has  a  given  probability  of  being  taken.  If 
either  activity  can  occur  prior  to  the  next  network  node,  the  form 
of  tne  probability  distribution  for  the  transit  time  to  next  node 
is  a  mixture  of  the  distributions  for  the  alternate  activity  times, 
with  weights  equal  to  the  probability  the  activity  is  taken.  This 
model  is  a  particular  Instance  of  a  semi-Markov  process,  a  type  of 
stochastic  process  frequently  observed  in  Industrial  operations. 

An  interesting  special  case  of  a  two- component  mixture  model  is 
one  in  which  the  components  (alternate  activities)  are  exponentially 
distributed.  The  form  of  the  c.d.f.  for  this  inter-node  duration  is 


F(t)  s  a(l  -  exp(-r1  t))  ♦  (1-a)(1  -  exp(-r2  t)), 


where  a  is  the  weight  associated  with  the  first  component,  and 
with  rate  parameters  r1  and  r2  for  the  1st  and  2nd  component  dist¬ 
ributions,  respectively.  The  p.d.f  for  this  mixture  model  is 


f(t)  s  a  rl  exp{-r1  t)  ♦  (1-a)  r2  exp(-r2  t). 


The  sum  of  n  such  "activities**  will  have  a  distribution  denoted 
by  g  (t),  for  the  p.d.f.,  and  by  G  (tl.  for  the  c.d.f.  of  time  t. 
n  n 

Using  the  convolution  theorem  of  Laplace  transforms,  the  transform 
of  g  (t)  can  be  written  as 
n 


g*(s)  r  (a  rl/(s  ♦  rl)  ♦  (1-a)  r2/(a  ♦  r2))  . 


T(i  facilitate  obtaining  an  inverse,  this  expression  Is  expanded 
in  i  power  series  of  terms  in 


i  n-i 

l/(s  ♦  rl)  /(s  ♦  r2) 
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The  mixed  products  in  this  series  must,  then,  be  expressed  in  a 
continued  fraction  expansion.  This  result  can  be  inverted  term  by 
term.  For  example,  for  n  =  2,  the  Laplace  transform  after  the 
indicated  operations  is 


2  2 

g*(s)  =  (a  r1)  /(s+r1)  +  2a(l-a)r 1r2/ (r 1-r2)/ (s^r2)  + 

2 

2  2 

((l-a)r2)  /(s+r2)  +  2a(  1-a)r  1r2/ (r2-r  1 )/ (s^-rl ) .  (20) 

The  inverse  transformation  is  obtained  by  inspection. 

2  2 

g  (t)  =  (a  r1)  t  exp(-rl  t)  +  ((l-a)r2)  t  exp(-r2  t) 

2 

+  2a(l-a)rlr2/(r1-r2) (exp(-r2  t)  -  exp(-rl  t))  .  (21) 

Integrating  g  (t)  produces  the  c.d.f.i 
2 


2  2 

G  (t)  =  1  -  a  (Url  t)  exp(-rl  t)  -  (1-a)  (Ur2  t)  exp(-r2  t) 

2 

-  2a(1.a)/(r1-r2)(r1  exp(-r2  t)  .  r2  exp(>r1  t)),  (22) 

Closed- form  expressions  for  G  (t)  for  larger  values  of  n  are  found 

n 

In  the  Implementing  computer  progrxm  in  Annex  D.  These  expressions 
are  used  to  calculate  the  Normal  approximation  error  for  the  c.d.f. 
Results  are  shown  in  Table  5  for  a  numerical  example  In  which  tne 
parameter  a  r  0.8,  and  the  mean  values  of  the  first  and  second  compon 
ents  are  in  the  ratio  of  0.05  to  1,0.  Rate  parameters  rl  and  r2  are 
adjusted  to  always  yield  a  mean  value  of  the  sum  equal  to  unity. 

The  last  practice  assures  that  the  same  points  are  evaluated  in  the 
domain  of  the  c.d.f.  regardless  of  the  value  of  n.  Also  note  that 
the  computer  program  (LP.INV)  uses  16  log- transformed  points  at  whlcn 
the  c.d.f.  error  is  evaluated — not  the  usual  20. 


TABLE  5 

ERROR  IN  THE  NORMAL  APPROXIMATION  OF  N  CONVOLUTIONS  OF  A 
TWO-COMPONENT  EXPONENTIAL  MIXTURE  PROD  DISTRIBUTION 


N 

Max  Abs  Error 

r"MS  Error  Off  points) 

2 

o'iSS 

0.  fsff 

3 

0.242 

0.166 

4 

0.213 

0.142 

5 

0.188 

0.124 

10 

0.117 

0.074 

15 

0.032 

0.052 

9 


The  maximum  absolute  errors  in  Table  5  are  nearly  3  times  the 
corresponding  errors  in  Table  3f  which  referred  to  convolutions 
of  a  single  exponential  component.  Further,  the  RMS  errors  in 
Table  5  are  about  4.0  to  4,5  times  the  corresponding  errors  in 
Table  3.  These  observations  indicate  that  the  sum  of  n  RVs  from 
an  exponential  mixture  distribution  may  converge  VERY  SLOWLY,  with 
increasing  n,  toward  Normality.  In  fact,  the  approximation  errors 
in  the  c.d.f.  of  the  sum  may  be  much  greater  than  comparable 
errors  in  sans  of  exponential  RVs,  (which  are  even  quite  large), 

A  frequently  used  rule  of  thumb  for  deciding  what  is  a  marginally 
large  sample  size  in  many  statistical  applications  is  that  n  >  30 
is  "large".  However,  for  the  30-fold  convolution  of  the  exponential 
mixture  c.d.f.,  one  finds  that  the  approximating  Normal  c.d.f.  has 
a  max  absolute  error  of  nearly  0.057  and  an  RMS  error  of  about  0.026. 
When  dealing  with  sums  of  RVs  from  a  semi-Markov  process,  consider¬ 
able  inaccuracy  can  be  encountered  in  taking  a  Normal  approximation. 
This  is  the  lesson  of  this  particular  example. 

12.  Monte-Carlo  Simulation 


As  is  shown  above,  closed- form  expressions  can  be  obtained  for 
convolutions  of  an  exponential  mixture  distribution  by  using  Laplace 
transform  methods.  However,  the  complexity  of  inverting  G»(s)  grows 

n 

rapidly  /ith  n.  In  this  case  in  particular,  alternatives  to  eval¬ 
uating  '  jrmulas  are  sought  for  calculating  points  of  G  (t)  for  large 

n 

values  of  n.  As  suggested  above  in  paragraph  7»  p*5f  Monte-Carlo 
is  a  useful  and  quite  general  technique.  For  example,  in  the  case 
of  the  exponential  mixture  model,  generation  of  one  RV  from  the 
mixture  distribution  involves:  U)  drawing  one  uniform  (0,1)  deviate, 
U;  (b)  drawing  a  RV  from  an  exponential  distribution  with  rate 
parameter  r1,  if  U  <  a;  or  (c)  otherwise,  drawing  a  RV  from  an 
exponential  distribution  having  rate  parameter  r2.  Tne  sum  of  n 
sucn  random  variables  is,  of  course,  the  RV  of  interest.  Run  time 
for  generating  an  estimate  of  G  (t)  by  simulation  is  actually  found 

n 

to  be  somewhat  less  than  that  indicated  by  equation  (13)  for  sums 
of  uniform  RVs,  due  to  a  different  choice  of  points  in  the  domain 
of  the  c.d.f.  at  which  the  distribution  is  evaluated.  The  particular 
numerical  example,  Introduced  in  paragraph  11,  is  used  to  compare 
a  Monte-Carlo  estimate  with  a  theoretical  estimate  and  with  a  Normal 
approximation  of  G  (t).  Results  for  two  values  of  n  are  displayed 
n 

in  Table  6.  Tne  theoretical  estimate  is  obtained  by  formula  eval¬ 
uation  for  n  s  3i  and  is  obtained  by  Bellman's  numerical  inversion 
method  (to  be  discussed),  for  n  s  IS*  For  this  problem  the  max  abs¬ 
olute  error  in  Bellman's  method  is  quite  small— typically  <  0.001 — 
making  tnis  a  good  theoretical  estimate.  Exponential  rate  parameters 
are  scaled  so  titat  the  mean  of  the  sum  Is  unity  for  both  values  of 
n.  The  effect  of  the  time  scaling  makes  the  variance  of  the  sum 
inversely  proportional  to  n.  Thus,  the  standard  deviation  of  the 
sum  is  1.U1SB  for  n=3i  and  is  0.6332  for  ns1S»  in  this  example. 

Note  that  tne  Normal  approximation  is  quite  poor  at  low  quantiles, 
even  for  n  as  large  as  15.  Also  note  that  the  Monte-Carlo  estimate 
is  quite  good  for  20,000  replications.  Tne  max  abs  error  is  nearly 
O.OOG,  and  the  RMS  error  is  about  0.003  for  one  random  number  stream. 

10 


TABLE  6 

SEVERAL  APPROXIMATIONS  OF  THE  C.D.F.  OF  THE  SUM  OF  N  RV’S 
FROM  A  TWO-COMPONENT  EXPONENTIAL  MIXTURE  PROB  DISTRIBUTION 


N  =  3  N  =  15 


Sum 

Value 

Theory 
Eval'  n 

Monte 

Carlo 

Norm 

Approx 

Theory 
Eval'  n 

Monte 

Carlo 

Norm 

Approx 

0.0053 

d.o'ooo' 

0.0000 

0.f4lf 

0.0600 

0.0000 

d.osffi 

0.0281 

0.0044 

0.0043 

0.2462 

0.0000 

0.0000 

0.0624 

0.0695 

0.0432 

0.0438 

0.2555 

0.0002 

0.0000 

0.0708 

0.1304 

0.1577 

0.1530 

0.2696 

0.0027 

0.0030 

0.0848 

0.2120 

0.3256 

0.3292 

0.2889 

0.0371 

0.0336 

0.1067 

0.3161 

0.4740 

0.4796 

0.3145 

0.1070 

0.1036 

0.1400 

0.i]450 

0.5668 

0.5696 

0.3475 

0.1956 

0.1938 

0.1904 

0.6024 

0.6218 

0.6259 

0.3894 

0.3109 

0.3088 

0.2650 

0.7930 

0.6646 

0.6676 

0.4419 

0.4473 

0.4478 

0.3718 

1.0239 

0.7073 

0.7123 

0.5067 

0.5948 

0.5929 

0.5150 

1.3057 

0.7519 

0.7542 

0.5855 

0.7360 

0.7330 

0.6854 

1.6552 

0.7981 

0.7996 

0.6782 

0.8545 

0.8571 

0.8496 

2.1013 

0.8451 

0.8473 

0.7817 

0.9370 

0.9381 

0.9590 

2.7003 

0.8918 

0.8949 

0.8851 

0.9817 

0.9811 

0.9964 

3.5859 

0.9367 

0.9394 

0.9661 

0.9974 

0.9974 

1.0000 

5.2401 

0.9771 

0.9764 

0.9986 

1.0000 

1.0000 

1.0000 

RMS  errors  in  the  Monte-Carlo  c.d.f.  estimate  seem  to  vary  inversely 
as  the  square  root  of  the  sample  size  (S)t  over  the  range  from  5  to  20 
thousand  replications,  and  do  not  vary  statistically  with  the  order 
(N)  of  the  convolution.  Typical  RMS  errors  for  this  example  over  thlJ. 
range  in  S  vary  from  0.002  to  O.OOU.  An  approximation  for  Monte-Carlo 
run  time  (sec)  on  the  Prime  9935  is 

T  5  S  N/^,  (23) 

where  S  is  given  in  thousands  of  replications,  and  with  2  le  N  le  20. 
Run  time— as  opposed  to  c.p.u.  time— is  dependent  on  the  number  of 
other  users  sharing  the  couputer  and  upon  the  nature  of  their  Jobs. 

The  value  of  T  given  here  is  representative  of  the  active  part  of  a 
work  day. 

13.  Bellman's  Metnod 


A  numerical  method  Is  given  In  Ref  la  by  Bellman  for  inverting 
Laplace  transforms.  Derivation  of  the  method  proceeds  from  the 
definition  of  the  Laplace  transform  of  an  analytic  function  F(t): 

F*{s)  =  lntegral(0,inf ):  exp(-st)  F(t)  dt.  (2*<) 

Tne  variable  of  integration  is  changed  to  x  via  the  transformation 

t(x)  s  In  (2/(x  ♦  1)).  (23) 
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Then,  the  real  variable  x  is  defined  on  (-1,1).  Now,  the  trans¬ 
form  variable,  s,  is  replaced  with  a  discrete  real  variable  c  k, 
with  c  constant  and  k  integer,  1  le  k  le  m.  That  is,  the  transform 
Is  evaluated  at  discrete,  evenly  spaced  points  on  the  real  line. 

The  variable  of  integration  is  also  discretized  Lt  m  points,  x(j), 

1  le  j  lem.  Tnus ,  the  integration  operation  is  replaced  by  a 
summation.  Gaussian  quadrature  is  chosen  as  the  mean  of  evaluat¬ 
ing  the  integral.  The  x(j)  are  chosen  as  the  points  of  the 
independent  variable  in  an  mth  order  quadrature.  Notationally,  let 

g(J)  =  F(t(x(j))),  j  =  1,  2,  ...,  ra.  (26) 

The  weight  function  [1]  for  gaussian  quadrature  is  denoted  by  w(j), 
for  1  le  j  le  m.  With  this  notation,  equation  (24)  becomes 

ck-1 

F»(ck)  =  Sum  over  J  (l,ra):  0.5  w(j)((x(j)+l)/2)  g(j),  (27) 

for  k  =  1,  2,  ,  m.  Tnis  equation  is  seen  to  be  a  matrix 

equation,  wliich  can  be  written  compactly  as 

F»  =  A  g  .  (28a) 


where  F”  and  g  are  m- component  column  vectors  and  where  a 

typical  element  a  of  the  A  matrix  is 
kj 

ck*1 

a  :  0.5  w(J)  ((x(J)4l)/2)  .  (28b) 

kJ 

Equation  (28)  is  solved  for  g  .  Then,  a  points  of  P(t)  are 

obtained  froo  (26),  with  associated  values  of  the  independent 
variable,  obtained  froo  (25)  •  For  the  best  accuracy  for 
the  c«d.f.  on  several  sample  problems  using  Bellman's  method, 

It  is  found  that  the  value  of  the  constant  c  should  be  unity 
and  that  the  problem  scale  parameters  should  be  adjusted  so 
that  the  mean  value  of  the  sum  (t)  Is  approximately  unity. 

(If  necessary,  rescaling  t  can  be  done  following  the  solution 
of  (28),  in  order  to  preserve  original  units  of  the  independent 
variable.)  It  Is  found  that  the  matrix  A  becomes  progressively 
closer  to  being  singular  as  m  Increases.  For  double- precision 
arithmetic  on  the  Prime  9955,  it  is  found  that  truncation  error 
Xislts  the  maximum  value  of  m  to  about  16.  However,  these  16 
points  of  Fd)  are  calculated  quite  rapidly  and  accurately.  For 
example,  the  HMS  error  To.**  3  convolutions  of  an  exponential  Is 
about  O.OOOOO^i,  and  the  RKS  error  for  3  convolutions  of  the  above 
^exponential  mixture  is  about  0*00027.  Thus,  when  the  Laplace 
tninsfom  of  a  distribution  of  interest  is  easily  and  accurately 
oalculated^  Del  I  man's  method  is  the  method  of  choice. 

Of  ^TheVei^ts  •  w(J),  and  the  points,  x(J),  for  mth  order  gaussian 
quadra t.ure  are  listed  in  Handbook  of  Mathematical  Functions, 


AMS  56  (1966),  on  page  916,  for  values  of  m  from  2  to  96. 

12 


14.  Convolutions  of  a  Two- Parameter  Wei bull  Distribution 


Analytic  methods  and/or  Bellman’s  inverse  Laplace  transform  are 
not  well  suited  to  obtain  convolutions  of  certain  types  of  prob¬ 
ability  distributions.  Examples  are:  (a)  a  distribution  with 
a  threshold  parameter  and  cii  upper  truncation  limit,  and  (b)  a 
distribution  whose  Laplace  transform  is  difficult  to  express  or  to 
evaluate  accurately.  One  probability  distribution  of  practical  [1] 
interest  which  suffers  from  the  last  difficulty  is  the  two- parameter 
Wei bull  function,  whose  c.d.f.  is  given  as 

b 

F(x)  =  1  -  exp(  -(x/a)  )  ,  0  le  x  <  inf,  (29) 

with  scale  parameter  a  and  shape  parameter  b. 

Altho  the  transform  can  be  expressed  as  an  error  function  of  s, 
the  result  is  difficult  to  evaluate  with  accuracy  sufficient  for 
inversion  via  Bellman’s  method.  Further,  closed-form  expressions 
for  the  n-fold  convolution  of  (29)  become  quite  complicated  for 
n  large.  The  closed- form  expression  for  n  =  2,  taken  from  Ref  1c, 
for  the  special  case  in  which  the  shape  parameter,  b,  s  2,  is 

2 

F  (t)  =  1  -  exp(-z  )  -  sqrt(pi/2)  z  (N(z)  -  N(-z)),  (30a) 

t 

where  N(z)  is  the  standard  Normal  integral  with  argument  z,  and 
with 

z  s  t/a  •  (30b) 

A  general  formula  which  approximates  the  n-foXd  convolution  of 
a  two- parameter  Weibull  distribution  was  derived  by  Leonard 
Johnson  (21.  The  LJ  approximation  is  an  Erlang  distribution  in 
the  argument  u,  where 

b 

u  s  (pt/a)  •  (31) 

Tne  parameter  p  is  chosen  so  that  the  mean  of  the  approximating 
distribution  matches  its  counterpart  in  the  convolution  distribution. 

p  5  gammaCn  ♦  l/b)/gamjia(l  ♦  1/b)/n!  ,  (32) 

with  complete  gamma  function  gammaiargument ). 

FT]  fheTwo^^rameter  Weibull  distribution  has  proved  to  be 
a  good  model  for  the  life  distribution  of  components  or 
systems  subject  to  fatigue  failure.  For  this  reason  it  is 
used  extensively  in  the  automotive  industry.  See  Ref  Ic. 

[2]  Johnson,  L.  GMR  Reliability  Manual,  GHR-302, 


General  Motors  Research  Labs,  i960. 


(33) 


Thus,  the  LJ  approximation  for  the  n-fold  convolution  c.d.f,  is 

i 

F  (t)  r  1  -  exp(-u)  Sum  over  t  (0,n-l)!  u  /i!  . 
t 

From  equations  (31,32,33)  it  is  seen  that  the  LJ  approximation 
is  exact  for  n  =  1 .  To  see  how  the  error  of  approximation  grows 
with  the  order  of  convolution,  consider  the  following  numerical 
example.  Let  the  scale  parameter,  a,  =  6  and  the  shape  parameter, 
b,  =  2.  To  compare  the  LJ  approximation  with  other  methods,  this 
problem  is  also  solved  using  these  other  methods;  (a)  discrete 
numerical  convolution,  using  1028  points  on  a  domain  that  comprises 
the  0th  to  the  99.97th  percentiles,  (b)  Normal  approximation,  and 
(c)  Monte-Carlo  simulation  with  a  sample  size  of  20,000  replications. 
The  RMS  error,  over  20  equi-spaced  points,  for  each  of  these  methods 
is  shown  in  Table  7.  The  error  for  the  discrete- numerical  (DN) 
convolution  Is  shown  for  n  =  2,  since  an  analytic  expression  exists 
as  a  check,  in  this  instance.  Since  this  error  is  relatively  quite 
small,  the  DN  solution  Is  used  to  evaluate  the  c.d.f.  errors  for  other 
values  of  n.  As  expected,  the  Normal  approximation  decreases  with  n. 
By  contrast,  the  LJ  approx  error  increases  with  n.  For  n  greater 
than  or  equal  to  7,  the  Normal  approximation  has  a  smaller  RMS  error 
than  the  LJ  approximation,  and  hence  is  preferred  to  LJ  there. 

RMS  errors  of  the  Monte-Carlo  (MC)  method  are  relatively  Independent 
of  convolution  order.  The  values  given  here  are  the  average  produced 
by  two  random  number  streams. 

TABLE  7 

RMS  ERRORS  IN  THE  C.D.F.  OF  THE  SUM  OF  N  IDENTICAL 

TWO-PARAMETER  WEIBULL  RV*S  PRODUCED  BY  SEVERAL  METHODS 


Convol *  n 

Method 

Of  Calculation 

Order  (N) 

DN 

LJ 

NA 

[ 

1 

0.0000 

o.oood"' 

6.o"i9U 

0.0013 

2 

o.ooou 

0.0028 

0.0145 

0.0024 

3 

0.0  (2) 

0.0036 

0.0100 

0.0022 

u 

0.0 

0.0001 

0.0081 

0.0026 

5 

0.0 

0.0005 

0.0067 

0.0026 

6 

0.0 

0.0009 

0.0060 

0.0028 

7 

0.0 

0.0050 

0.0050 

0.0036 

Li)  Average  value  of  the  error  over  two  random  number  streams. 

[2]  Value  of  the  discrete  numerical  error  Is  not  evaluated  for 

n  >  2,  but  Is  considered  relatively  small  versus  other  errors. 

lb.  Run  Time  Comparisons 


Wnereas,  the  discrete  numerical  method  is  quite  accurate,  and  is 
rensonubly  fast  for  n  =  2,  run  time  for  this  method  Increases  as  a 
power  function  of  n  -  2,  with  a  power  of  about  1.25.  Thus,  for  a 
constant  density  of  \02U  points  on  the  domain  of  the  convolution 
c.d.f.,  an  approximate  run  time  (sec)  is  given  by 


T  =  40  (n.2) 


1.25 


(34) 


For  n  =  2  only  20  points  of  the  numerical  convolution  are  evaluated. 

For  this  reason  run  tine  is  a  fraction  of  a  second  for  n  r  2,  whereas 
for  larger  values  of  n,  the  maximum  number  of  points  (1024)  are  calc¬ 
ulated  for  each  of  the  convolutions  except  the  last  (nth).  It  is 
emphasized  that  run  time  for  any  method  strongly  depends  upon  the 
background  activity  of  the  (time-shared)  computer.  Equation  (34)  gives 
nearly  the  maximum  time  experienced.  Minimum  run  times  are  approx¬ 
imately  half  of  maximum.  The  form  of  equation  (34)  suggests  that 
computational  overhead,  e.g.  paging,  increases  faster  than  n  does. 

When  n  is  3f  the  Monte-Carlo  run  time  for  20,000  replications  is 
about  the  same  as  the  run  time  for  the  DN  method.  However,  for 
higher-order  convolutions,  MC  is  faster.  For  example,  for  n  r  4, 

DN  requires  50jt  more  time  to  execute  than  MC.  For  nr  5,  DN 
requires  75X  more  time  to  run  than  MC,  i.e.,  the  ratio  of  run  times  is 
about  1.75.  For  n  r  8,  this  ratio  is  2.6.  Thus,  if  one  is  satisfied 
with  an  RMS  error  less  than  0.35tf  Monte-Carlo  would  be  the  preferred 
of  these  two  methods,  for  n  >  3.  Considering  the  errors  of  the  LJ 
and  NA  methods,  these  are  not  very  attractive  unless  execution  time 
is  a  major  consideration.  If  minimum  run  time  is  a  primary  consider¬ 
ation  for  this  type  of  problem,  a  hybrid  method  might  be  used  in  which 
DN  is  used  for  n  <  4,  LJ  used  for  4  le  n  <  7,  and  NA  used  for  n  ge  7, 
The  computer  source  program  (INT.TEST)  used  in  making  the  comparisons 
in  Table  7  is  found  in  Annex  E. 


16.  Fourier  Transform  Method 


As  noted  above  (p.  3»  5)»  the  product  of  an  integral  transform 

of  each  of  two  functions  corresponds  to  the  transform  of  the  convol¬ 
ution  of  the  functions.  This  theorem  has  already  been  exploited  In 
connection  with  the  Laplace  transform.  Tnls  paragraph  is  concerned 
with  an  application  of  this  theorem  using  the  Fourier  transform. 

An  important  and  practical  Fourier  transform  method  uses  an  algoritnm 
for  calculating  the  finite  Fourier  transform  (or  its  inverse)  due  to 
Cooley  and  Tukey,  and  called  the  fast  Fourier  transform  or  FFT  (1). 

The  speedy  execution  of  the  FFT  makes  practical  the  following  method. 
Two  density  functions  are  each  evaluated  at  a  particular  number  of 
equi-spaced  points  on  their  domains.  Tnese  data  vectors  are  input  to 
the  FFT,  which  yields  the  complex- valued  transforms.  These  transforms 
are  multiplied  (observing  the  rules  of  complex  arithmetic)  to  obtain 
the  transform  of  the  convolution  density.  Finally,  the  inverse  FFT  is 
performed  on  this  function  to  yield  the  required  density.  In  the 
computer  program  for  performing  these  operations,  found  In  Annex  E,  a 
function  f(x)  is  represented  in  complex  form  by  a  set  of  n  points  In 
which  there  are  n/2  real  components  and  n/2  imaginary  components.  Note 
that  n  must  be  an  Integer  power  of  2  for  this  purpose.  These  real  and 
complex  components  are  stored  in  adjacent  storage  locations  in  tne 
n-elcment  vector.  Of  course,  the  densities  being  convolved  have  only 
real  components,  so  that  all  imaginary  components  of  f(x)  are  assigned 
0  value.  Since  the  transform  occurs  in  place,  the  transform  of  f(x). 


[1)  Bloomfield,  P.  Fourier  Analysis  of  Time  Series:  An  Introduction, 
John  Wiley,  New  York,  NY,  c.  1976. 
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denoted  as  f»(w),  is  also  stored  in  the  n-element  vector  with  real  and 
Imaginary  components  of  the  transform  also  located  in  adjacent  posi¬ 
tions,  In  general,  the  Imaginary  components  of  f*(w)  are  non-zero. 

The  n/2  real  frequency  components  of  f»(w)  are  denoted  by  f»(w(k)) 
with  k  odd,  and  the  n/2  imaginary  frequency  components  are  located 
In  elements  of  the  vector  f*(k)  with  k  even,  (k  =  1,  2,  n). 

In  this  formulation  the  transform  and  its  inverse  are  duals  related 
by  equations  (35)  ard  (36)t 

f»(w(k))  =  Sum  over  j  (1,n)t  exp  (-iw(k)(j-1) )  f(x(j))/n,  (35) 
where  w(k)  is  the  kth  complex  frequency  with 

w(k)  s  2  pi  (k-1)/n  ,  k  r  1,  2,  n, 

and  where  I  is  the  pure  imaginary,  sqrt(-l). 

Then, 


f(x(j))  s  Sum  over  k  (l,n)t  exp  (iw(k)(j-l))  f«(w(k)),  (36) 

Because  of  the  dual  nature  of  f(x)  and  f*(w),  the  same  routine 
that  produces  a  transform  can  obtain  an  inverse  transform  merely  by 
specif ing  which  type  of  operation  is  wanted  via  "sign**  s  -1  for  the 
Fourier  transform,  and  by  sign  *  1  for  an  inverse  transform.  The 
computer  code  for  this  algorithm  is  found  in  Annex  F, 

17,  A  series  of  numerical  tests  were  performed  for  accuracy  and 
run  time  using  the  Fourier  transform  method.  These  are  compared 
with  Honte-Carlo  tests  using  the  same  test  functions*  Probability 
densities  used  as  test  functions  have  the  standardized  Erlang  and 
standardized  Wei  bull  forms.  In  both  instances  the  scale  parameter 
is  unity,  and  the  function  is  characterised  by  Just  a  shape  para¬ 
meter,  In  the  first  numerical  example  with  n  Erlang  densities 
being  convolved,  n-1  of  these  have  been  assigned  a  shape  parameter 
of  2  and  one  is  assigned  shape  parameter  3*  RHS  errors  are  shown  in 
Table  8,  for  selected  values  of  n,  for  both  the  Fourier  transform  (FFT) 
method  and  for  a  Monte-Carlo  simulation  with  20,000  replications. 

Tne  RMS  error  is  obtained  over  16  equi-spaoed  points  on  the  domain. 

The  number  of  points  (equivalently,  real  Fourier  frequencies)  used 
to  represent  the  densities  is  also  a  parameter  in  these  tests, 

TABLE  8 

RMS  ERRORS  IN  THE  C.D.F*  OF  THE  SUM  OF  N  RV'S  FROM 
7^  ERLANG  DISTRIBUTIONS  VIA  FFT  AND  HONTE-CARLO  METHODS 


Con  V  Jl  *n  ypT  ViTh  V  Vea  f  Tr^quenci  es  "Monte 


Order  (M) 

1024 

20«8 

4096 

Carlo  (201:  reps) 

2 

o.offoff 

0,0004 

*0.15352 

0.332? 

3 

0,0016 

0.0008 

0.0004 

0.0019 

4 

0,0023 

0.0012 

0.0006 

0.0007 

5 

0,0031 

0.0016 

0.0009 

0.0010 

10 

0,0072 

0.0036 

0.0018 

0.0008 

20 

0,0136 

0.0032 

0.0041 

0,0011 
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Using  the  FFT  method  with  2048  real  frequencies,  the  run  time  for 
any  value  of  n  varies  from  about  9  to  18  seconds.  Run  time  for  the 
FFT  method  seems  to  be  dominated  by  the  time  to  obtain  the  Fourier 
transforms  and  to  obtain  the  inverse.  Because  relatively  little  time 
Is  spent  in  multiplying  transforms,  run  time  is  essentially  independ¬ 
ent  of  the  convolution  order  for  the  values  shown.  By  contrast,  it  Is 
seen  that  Monte-Carlo  run  time  (T)  increases  nearly  linearly  with  N: 

T  =  8(N  -  2)  ♦  12.  (37) 

(This  approximation  is  quite  similar  to  that  given  in  equation  (13) 
for  evaluating  the  c.d.f,  of  the  sum  of  N  uniform  random  variables 
at  20  points  via  Monte-Carlo.)  Run  time  for  the  FFT  method  does 
increase  in  a  proportional  manner  with  #  of  real  Fourier  frequencies. 
For  example,  when  1024  real  frequencies  are  used  run  time  is  about 
5  seconds.  This  time  increases  to  9  seconds  for  2048  real  frequencies 
and  to  about  20  seconds  for  4096  real  frequencies.  Thus,  in  terms 
of  run  time,  calculating  the  c.d.f.  of  the  sum  of  3  Erlang  RVs  is 
nearly  the  same  using  either  Monte-Carlo,  with  20,000  replications, 
or  the  FFT  method,  using  4096  real  frequencies.  It  is  noted  that  for 
a  high-order  convolution  Integral— say,  >  10— a  very  large  number 
of  Fourier  frequencies  are  required  to  make  the  accuracy  of  the  FFT 
method  competitive  with  Monte-Carlo.  This  point  is  illustrated  by 
the  results  in  Table  8.  It  is  also  demonstrated  by  another  numerical 
example.  Consider  the  case  in  which  all  the  distributions  being 
convolved  are  standardized  exponential.  The  RMS  errors  for  the  FFT 
and  Monte-Carlo  methods  for  this  case  are  shown  in  Table  9.  Note 
that  these  results  are  substantially  the  same  as  those  in  Table  3. 

TABLE  9 

RMS  ERRORS  IN  THE  C.D.F.  OF  THE  SUM  OF  N  IDENTICAL 

EXPONENTIAL  RV*S  USING  FH  AND  HONTE-CARLO  METHODS 


ConwoV^n 

Numerical 

Method  • 

Order  (N) 

Fourier  Transform 

Monte  Carlo 

. i . 

. ~576SSi^~'~ 

6.66i6 

3 

0.00045 

0.0018 

4 

0.00064 

0.0019 

5 

0.00032 

0.0022 

10 

0.00178 

0.0010 

•  The  FFT  method  Implemented  here  has  4096  real  Fourier  frequencies 

(8192  element  array).  Simulation  sample  is  20,000  replications. 

Monte-Carlo  results  shown  are  averages  for  two  random  streams. 

18.  A  t  ilrd  numerical  example  was  used  to  test  tne  accuracy  of  the 
Fourier  . 'ansform  method.  In  this  case  a  standardized  o^eibuU 
density  w  th  shape  parameter  z  2  is  convolved  n  times  to  yield  the 
p.d.f.  for  ihe  sura  of  n  such  Welbull  RVs.  For  the  particular  case 
in  which  n  is  2,  the  numerical  error  in  the  c.d.f.  is  found  by 
comparing  the  exact  result  frora  equation  (30)  with  tne  FFT  approx¬ 
imation.  The  RMS  error  for  this  case  Is  0.00089#  about  three 
times  that  for  the  previous  two  examples.  Tnus,  the  numerical  error 
of  the  Fourier  transform  method  Is  rather  sensitive  to  the  form  of 


17 


the  distributions  being  convolved.  Simpson's  rule  is  used  to  calc¬ 
ulate  mean  and  SD  from  a  numerical  c.d.f.  for  the  sum  of  two  Wei bull 
RVs  from  the  distribution  with  shape  parameter  s  2.  For  the  case  of 
4096  real  Fourier  frequencies,  the  error  in  the  FFT  mean  is  0.186)1. 

By  contrast,  the  error  In  the  mean  value  using  the  discrete  numerical 
(DN)  convolution  with  1024  points  Is  0.09% •  A  comparable  relationship 
exists  in  the  RMS  error  In  the  c.d.f.  for  the  DN  versus  the  FFT  method. 
RMS  error  in  the  convolution  c.d.f.  for  this  example  Is  0.00036,  for 
DN,  versus  0.00089  for  FFT.  By  either  measure  of  error,  the  DN 
method,  applied  on  a  set  of  1024  points,  Incurrs  less  than  half  the 
error  of  the  FFT  method,  applied  on  a  set  of  4096  points.  If 
accuracy  of  results  were  tne  sole  criterion,  discrete  numerical 
convolution  would  certainly  be  preferred  to  the  FFT  method.  However, 
for  high-order  convolutions,  DH  Is  computationally  expensive  relative 
to  FFT.  For  example,  5  convolutions  of  a  Wei bull  distribution  using 
DN  with  this  degree  of  discretization  takes  about  160  sec.  (equa¬ 
tion  (34)),  In  a  comparable  run  environment,  FFT  with  4096  real 
Fourier  frequencies  requires  about  50  sec  for  the  same  problem. 

Thus,  the  FFT  method  executes  this  problem  In  one  third  the  time 
required  by  the  DN  method,  given  the  specified  density  of  points. 

It  is  noted  that  the  maximum  number  of  real  frequencies  (4096)  used 
with  FFT  in  the  above  examples  Is  the  maximum  permitted  on  our  Prime 
computer.  The  computer  system  limit  on  the  number  of  double- precision 
words  allocated  to  a  vector  Is  less  than  16,384.  If  the  number  of 
real  frequencies  were  doubled,  to  8192,  the  dynamic  storage  required 
for  both  real  and  iotaglnary  frequency  components  would  be  16,384. 

19.  Summary  and  Conclusions 


This  report  has  surveyed  several  methods  for  calculating  probability 
distributions  of  sums  of  independent  random  variables.  Formulas  for 
tne  c.d.f.  of  the  sum  have  been  derived  for  several  cases.  These 
cases  include  n  random  variables  fromi  (a)  a  standard  uniform  dist¬ 
ribution,  (b)  uniquely  different  uniform  distributions,  (c)  an  Erlang 
distribution.  Including  the  exponential  as  a  special  case,  (d)  dif¬ 
ferent  exponential  distributions,  (e)  two-component  exponential  mix¬ 
tures,  and  (f)  a  Welbull(2)  distribution  (two  RVs  only).  The  closed- 
form  solutions  were  used  to  evaluate  the  accuracy  of  various  numer¬ 
ical  methods,  including  approximations. 

20.  The  sum  of  n  RVs  from  some  distributions  have  a  c.d.f.  which 
rapidly  approaches  Normality  with  increasing  n.  Examples  of  this 
sort  are  the  uniform  distribution  ani  distributions  which  appear 
Normal,  such  as  gamma  with  large  shape  parameter.  However,  other 
distributional  forms  exhibit  relatively  slow  convergence.  These 
include  exponential  and  exponential  mixture  distributions.  The  last 
is  particularly  slow  in  converging  toward  Normality.  For  this  case, 
the  sum  of  15  RVs  has  a  c.d.f.  whose  Normal  approximation  nas  a  max 
aosolute  error  of  more  than  0.08,  which  is  intolerably  large  for  most 
purposes.  Generally,  the  c.d.f.  errors  of  a  Normal  approximation 
are  larger,  for  a  given  n,  if  the  scale  (or  rate-)  parameters  of  the 
component  distributions  exhibit  a  large  range  than  if  all  distributions 
are  identical. 

21.  It  is  difficult  to  make  unqualified  statements  concerning  the 
superiority  of  any  one  of  the  nuaerlcal  methods.  This  situa¬ 
tion  is  due  in  part  to  the  diversity  of  user  requirements  for  speed 
and  accuracy  and,  in  part,  to  the  fact  that  some  methods  are  parti- 
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cularly  suited  to  just  some  classes  of  distributions.  For  example, 
Bellman’s  method  requires  that  the  Laplace  transform  of  the  probab¬ 
ility  density  be  easily  and  accurately  calculated.  The  class  of  gamma 
distributions  and  of  mixtures  of  gamma  distributions  are,  therefore, 
well  suited  to  this  method.  Because  of  superior  run  time  and  accur¬ 
acy  [1],  Bellman's  method  is  the  method  of  choice  for  distributions 
in  the  gamma  class,  when  certain  conditions  are  met.  These  are: 

(a)  the  user  must  be  satisfied  with  a  sixteen-point  characterization 
of  the  p.d.f.  and  of  the  c.d.f.,  and  (b)  available  macnine  arithmetic 
can  operate  on  a  floating-point  word  with  48  bit  mantissa  and  7  bit 
(or  more)  exponent.  (These  requirements  are  met  on  the  Prime  9985 
minicomputer  with  double-precision  arithmetic.)  The  last  feature 
is  necessary  to  avoid  truncation  error,  which  is  critical  to  Bellman’s 
method.  In  those  instances  where  Bellman's  method  is  inapplicable, 
discrete  numerical  convolution  of  the  component  distributions  offers 
the  greatest  potential  for  accuracy,  at  a  cost  of  run  time.  Where 
run  time  is  an  important  consideration  as  well  as  accuracy,  use  of  tne 
Fourier  transform  method  with  the  FFT  algorithm  is  attractive,  pro¬ 
viding  the  order  of  convolution  does  not  exceed  about  ten.  (This 
statement  presumes  that  the  max  vector  dimension  <  16,384.)  Another 
advantage  of  the  Fourier  transform  method  is  that  it  is  quite  flexible 
with  regard  to  distributional  forms  that  can  be  handled.  Of  course, 
the  Normal  approximation  is  preferred  in  those  instances  where  the 
form  of  the  component  distributions  assures  rapid  convergence  toward 
Normality.  For  distributions  on  a  bounded  domain,  such  as  the  uniform, 
relatively  small  RMS  errors  in  the  c.d.f.  by  Normal  approximation  are 
incurred  when  the  order  of  convolution  Is  5  or  more.  In  cases  where 
accuracy  is  not  too  stringent — say,  an  RMS  error  of  0.002 — Monte-Carlo 
simulation  [2]  is  the  most  flexible  and  resonably  efficient  method 
studied.  A  somewhat  surprising  finding  is  that  Monte-Carlo  is  pre¬ 
ferred,  in  many  cases,  to  discrete  numerical  convolution  when  the 
tolerable  RMS  error  is  about  0,2%  and  when  the  number  (n)  of  random 
variables  in  the  sum  is  three  or  more.  Monte-Carlo  run  time  increases 
linearly  with  n,  but  the  rate  of  increase  is  not  as  great  as  that  for 
discrete  numerical  convolution.  In  comparing  Monte-Carlo  with  FFT,  it 
is  noted  that  the  RMS  error  for  Monte-Carlo  does  not  increase  with  n, 
as  the  FFT  error  does.  When  limited  by  computer  storage  to  4095  real 
frequency  components,  the  FFT  method  becomes  less  accurate  than 
Monte-Carlo  for  n  greater  than  about  ten.  Also,  time  to  code  a  given 
application  for  a  Monte-Carlo  simulation  is  generally  the  least  of  any 
method. 


flT  ^Numer  leaf ’error  in  the  c.d.f.  of  convolutions  of  tne  gamma  family 
have  error’s  via  Bellman's  metnod  of  the  order  of  10**-5. 

[2]  A  Monte-Car’lo  sample  of  20,000  replications  was  used  for  nearly 
all  numerical  tests.  This  sample  size  is  a  practical  value  in 
view  of  these  facts:  (a)  run  time  is  propoi’tional  to  sample  size, 
and  (b)  RMS  error  Is  Inversely  proportional  to  square  root  of  tne 
sample.  Halving  the  RMS  error  would  increase  run  time  by  a  factor 
of  4.  For  an  RMS  error  in  tne  c.d.f.  of  much  less  than  0.21,  tne 
required  Monte-Carlo  run  time  would  make  this  method  non¬ 
competitive  with  others. 
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DISTRIBUTION 


Copies 


1  HQDA  DAMO-ZA 

WASH  DC  20310 

1  HQDA  DALO-SMZ 

WASH  DC  20310 


COMMANDER 

USAMC 

5001  EISENHOWER  AVE. 

ALEXANDRIA,  VA  22333-0001 
1  ATTN:  AMCRE-IP 

1  AMCPA-S 

1  AMCDP 


DIRECTOR,  US  AMSAA 
ABERDEEN  PG,  MD  21005-5066 
1  ATTN:  AMXSY-DL 

1  AMXSY-R 

1  AMXSY-MP 


COMMANDER  USA 

C  ,MMUNICATIONS  AND  ELECTRONICS 
COMMAND 

FT  MONMOUTH,  NJ  07703-5304 
1  ATTN;  Al^EL-PL-SA 


COMMANDER  CECOM  (R&D) 

FT  MONMOUTH,  NJ  07703-5304 
1  ATTN:  AMSEL-SAD 


COMiiANDER  USAMICOM 
REDSTONE  ARSENAL, 

AL  35BO9-5O0O 

1  ATTN:  AMSMI-DS 

COMiMANDER  U3ATAC0M 
WARREN,  MI  48090 

1  ATTN:  AMSTA-V 

OFFICE  OF  PROJECT  MGR 
CANNON  ARTY  WPNS, DOVER 
NJ  07301-5001 

1  ATTN:  A.CPM-CAWS 

COMMANDER,  US  ARMY  LOGISTICS  CENTER 
FORT  LEE,  VA  23801 
1  ATTN;  ATCL-S 
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1  COMMAHDER 

DEFENSE  LOGISTICS  STUDIES 
INFORMATION  EXCHANGE 
FORT  LEE,  VA  23801 

C0:“1MANDSR 

USA  LOGISTICS  EVAL  AGENCY 
NEW  CUMBERLAND  ARMY  DEPOT 
NEW  CUMBERLAND,  PA  17070 
1  ATTN:  DAL-LEM 

COMMANDER 
US  M3SA 

LEXINGTON,  KY  110511-5101 
1  ATTN:  AMXMD-ER 

DIRECTOR,  US  ARMY 
INVENTORY  RESEARCH  OFFICE 
ROOM  800,  CUSTOM  HOUSE 
2  ND  &  CHESNUT  STREETS 
PHILADELPHIA,  PA  19106 
1  ATTN:  AMXMf'-IRO 

COMMANDER  USATECOM 
ABERDEEN  PG,  MD  21005-5055 
1  ATTN:  AMSTE-SY 

12  DEFENCE  TECHNICAL  INFORMATION  CENTER 
CAMERON  STATION 
ALEXANDRIA,  VA  22314 


COMMANDER  US  ARDEC  (D) 
DOVER,  NJ  07801-5001 


1  ATTN:  SMCAR-  -LC  (D) 

1  -SC  (D) 

1  -SE  (D) 

1  -RAA  (D) 

1  -MSI  (D) 


COMMANDER  US  AMCCOM  (R) 
ROCK  IS,  IL  61299-6000 
1  ATTN:  AMSMC- 
1 
1 
1 
1 
7 
1 


-AS  (R) 

-IR  (H) 

-QA  (R) 

-MA  (R) 

-OP  (H) 

-SA  (R) 
-IMP-L  (R) 
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DIRECTOR,  AMCCOM 
AMMO  CENTER 
SAVANNA,  IL  61074 

1  ATTN:  SARAC-DO 

COI^ANDER 

WATERVLIET  ARSENAL 
WATERVLIET,  NY  12189-5000 
1  ATTN:  AMXMC-LCB-TL 

COMMANDER 

CHEMICAL  R  AND  D  CENTER 
ABERDEEN  PROVING  GROUND 
(EDGEWOOD  AREA),  MD  21010-5423 
1  ATTN:  AMSMC-CLJ-IA  (A) 

DIRECTOR  US  AMETA 
ROCK  IS,  IL  61299-6000 
1  ATTN:  AMXOM-QA 

DIRECTOR 

NAVAL  POSTGRADUATE  SCHOOL 
MONTEREY,  CA  93940 
1  ATTN:  DEPT  OF  OPERATIONS  ANAL. 

1  DIRECTOR 

ADVANCED  RESEARCH  PROJECTS  AGENCY 
1400  WILSON  BLVD 
ARLINGTON,  VA  22209 

DIRECTOR 
USA  TRASANA 

WHITE  SANDS  MISSILE  RANGE 
WHITE  SANDS,  NM  88002-5502 
1  ATTN:  ATAA-SL 

COMMANDER 

USA  COMBINED  ARMS  COMBAT 
DEVELOPMENT  ACTIVITY 
FT  LEAVENWORTH,  KS  66207 
1  ATTN:  ATZL-CAM-M 


COMPUTER  SOURCE  PROGRAMS 


Source  programs  listed  in  Annexes  A-F  are  written  in  SIMSCRIPT  2.5 
for  the  PRIME  minicomputer.  However,  the  source  code  does  not  employ 
features  peculiar  to  this  computer.  Each  Annex  contains  a  MAIN  or  ex¬ 
ecutive  program  and  several  routines  and  functions.  At  the  beginning 
of  each  program  listing  are  found  a  functional  description  and  an  I/O 
list.  All  utility  functions  and  routines  are  included  among  these 
listings.  Inputs  to  MAIN  programs  are  read  interactively,  with  prompt¬ 
ing  messages  displyed  at  the  terminal.  No  external  files  are  used. 
Since  output  is  lengthy,  it  is  necessary  to  set  up  a  COMO  file  to 
display  all  of  it  and  to  obtain  a  permanent  copy.  The  functions  of 
the  MAIN  programs  are  summarized  here: 

RUN.NFOLD.U,  in  Annex  A,  obtains  the  probability  density  and  cum¬ 
ulative  probability  distribution  of  the  n-fold  convolution  of  a  stand¬ 
ard  uniform  distribution.  Approximations  to  the  p.d.f.  and  c.d.f.  of 
convolution  based  on  a  Normal  probability  function  are  calculated  and 
printed  for  comparison, 

RUN.NFOLD.GU ,  in  Annex  B,  calculates  and  prints  the  p.d.f.  and 
c.d.f.  of  the  sum  of  a  set  of  n  uniform  random  variables  drawn  from 
distributions  having  a  common  threshold  parameter  but  with  different 
domains.  Normal  probability  approximations  to  the  p.d.f,  and  c.d.f, 
are  calculated  and  printed  for  comparison  with  exact  results.  The 
maximum  absolute  error  and  the  RMS  error  are  calculated  for  the  Normal 
approximation  to  the  c.d.f.  Optionally,  a  Monte-Carlo  simulation  can 
be  performed  and  error  statistics  calculated  and  printed. 

RUN.NFOLD.E,  in  Annex  C,  obtains  the  p.d.f,  and  c.d.f,  of  the 
sum  of  n  exponential  random  variables.  Two  options  are  available; 

(a)  all  exponential  random  variables  are  from  the  same  dlstrlbfitlon, 
and  (b)  each  exponential  RV  is  from  a  uniquely  different  distribution. 
The  exact  c.d.f.  is  compared  with  a  Normal  approximation  on  a  finite 
point  set.  Max  abs  and  RMS  errors  are  calculated  and  printed. 

LP.INV,  In  Annex  D,  obtains  the  p.d.f.  and  the  c.d.f.  of  the 
sum  of  n  random  variables  from  Erlang  distributions  and  exponential 
mixture  distributions.  A  numerical  method  based  upon  the  Laplace 
transform  is  used  to  obtain  approximate  results.  This  method  involves 
calculating  the  Inverse  transform  via  Bellman’s  method.  A  closed-form 
solution  to  the  problem  Is  used  to  calculate  the  error  in  Bellman's 
method  and  the  errors  of  a  Normal  approximation  and  of  a  Monte-Carlo 
estimate  of  the  c.d.f. 

INT.TEST,  in  Annex  E,  tests  a  variety  of  methods  for*  obtaining 
convolution  Integrals  of  a  two- parameter  Weibull  dlstrl bution.  Tne 
methods  being  compared  are:  (a)  evaluation  of  an  analytic  expression, 

(b)  Leonard  Johnson’s  approximation  based  on  the  Erlang  distribution, 

(c)  discrete  numerical  convolution,  and  (d)  Monte-Carlo  simulation. 

The  max  absolute  error  and  the  RMS  error,  over  a  finite  set  of  points, 
are  calculated  and  printed  for  each  namerical  approximation. 

TEST.CONVOLV,  in  Annex  F,  obtains  convolutions  of  either  standard¬ 
ized  Erlang  or  Weibull  distributions  using  a  numerical  method  based  on 
the  finite  Fourier  transform.  Comparisons  with  exact  results  and, 
optionally,  with  Monte-Carlo  estimates  arc  also  given. 


ANNEX  A 


SIMSCRIPT  SOURCE  PROGRAM:  RUN.NFOLD.U 


1  PREAMBLE  "RUN.NFOLD.U 

2  NORMALLY  MODE  IS  REAL 

3  DEFINE  SNORM  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

4  DEFINE  ERRFX  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

5  END  "PREAMBLE 

1  MAIN  "RUN.NFOLD.U 

2  •  • 

3  "Driver  program  to  obtain  the  probability  density  and  cum  probability 

4  "of  the  N-fold  convolution  of  a  standard,  uniform  dlst.  PDFs  and  CDFs 

5  "of  the  Normal  dlst  having  same  mean  and  SD  Is  printed  for  comparison. 

6  " 

7  DEFINE  I,N  AS  INTEGER  VARIABLES 
6  DEFINE  ANSWER  AS  A  TEXT  VARIABLE 
9  'LO'SKIP  1  LINE 
10  PRINT  4  LINES  THUS 

This  program  calculates  and  prints  the  p.d.f.  and  c.d.f  of  the  N-fold 
convolution  of  a  standard,  (0-1)  uniform  probability  distribution.  User 
Inputs  are  the  Integer  N  and  the  upper  probability  limit  (PMAX)  to  terminate. 

15  PRINT  1  LINE  THUS 

INPUT  THE  VALUE  OF  N. 

17  READ  N 

18  PRINT  1  LINE  THUS 

INPUT  THE  (MAX)  VALUE  OF  THE  CDF  TO  TERMINATE  CALCULATIONS. 

20  READ  PMAX 

21  LET  AVG5N/2.O 

22  LET  VARsN/12.0 

23  LET  C0ND=1.0/SQRT.F(2.0»PI.C»VAR)  "FOR  NORMAL  DENSITY  COEF 

24  LET  STDV=SQRT.F(VAR) 

25  LET  LIMsAVG  ♦  3.0»STDV 

26  LET  LIM=MIN.F(REAL.F(N),  LIM) 

27  LET  DELTsLIM/20.0 

28  LET  LINES. V=9999 

29  SKIP  2  LINES 

30  PRINT  6  LINES  WITH  N 

31  THUS 

PROBABILITY  DISTRIBUTION  OF  A  ••-FOLD  CONVOLUTION  OF  A  STD  UNIFORM  DISTRIBUTION 


Indep 

Variable 

N-foXd  Convolution  Norail 

p.d.f.  c.d.f.  p.d.f. 

Prob  Distrlb 
c.d.f. 

”“37* 

LET  MAEsO.O 

39 

LET  RMSrO.O 

40 

FOR  1=1  TO  20  DO 

41 

LET  T=I»DELT 

42 

CALL  NFOLO.U  (N,T)  TIELDIHC 

PDF, CDF 

43 

LET  ARS:(T-AVG)/STOV 

44 

LET  NPDF=C0MD»EXP.F(-0.5»AB0«»2) 

45 

LET  NCOF^SSORHCARG) 

46 

LET  DIFF=CDF-NCDF 

Difference 
c.d.f . 
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47  LET  MAE=MAX.F(MAE,ABS.F(DIFF)) 

48  ADD  DIFF»»2  TO  RMS 

49  PRINT  1  LINE  WITH  T,  PDF,  CDF, 

50  THUS 

»,»»»»»»  tt.ttttttttttii 

52  IF  CDF  GE  PMAX 

53  GO  TO  LI 

54  OTHERWISE 

55  LOOP  ’’OVER  I 

56  'Ll ’PRINT  2  LINES  THUS 


NPDF,  NCDF,  DIFF 

».»•»»»»  ». 


59  LET  RMS=SQRT.F(0.05»RMS) 

60  PRINT  2  LINES  WITH  MAE, RMS 

61  THUS 

Max  ab3  error  in  Normal  approximation  of  c.d.f .  _  »,»»»»»» 

RMS  error  in  Normal  approximation  of  c.d.f.  _ 

04  PRINT  1  LINE  THUS 

DO  YOU  HAVE  OTHER  VALUES  OF  N?  (YES  OR  NO). 

60  READ  ANSWER 

67  IF  SUBSTR.FC ANSWER, 1,1)  =  ”Y” 

68  GO  TO  LO 

69  OTHERWISE 

70  STOP 

71  END  "MAIN 

1  ROUTINE  NFOLD.U  GIVEN  N,  T  YIELDING  PDF,  CDF 

2  ” 

3  ’’Routine  calculates  the  probability  density  function  (PDF)  and  the  cum* 

4  ’’ulative  distribution  function  (CDF)  of  the  N*foid  convolution  of  a 

5  ’’standard  uniform  (0,1)  probability  dist.  Real-valued  argument  is  T. 

6  ” 

7  ’’With  the  following  notation  for  the  CDF  argument  ti  F(n,t),  with 
3  ’’the  combination  of  n  things  taken  i  at  a  tine  denoted  as  C(n,i), 

9  ’’and  with  the  unit  step  function  at  x  denoted  by  u(t-x), 

10  ” 

11  ”  i  n 

12  ”  F(n,t)  5  Sum  (i=0  to  n):  (-1)  C(n,i)  u(t-i)  (t-i)  /n! 

13  *• 

14  DEFINE  I,N  AS  INTEGER  VARIABLES 

15  IF  T  LE  0.0 

16  LET  PDrsO.O 

17  LET  CDFsO.O 

1 8  RETURN 

1.)  OTHERWISE 

30  IF  T  GE  REAL.F(N) 

21  LET  PDFrO.O 

22  LET  CDF 5 1.0 

23  RETURN 

24  OTHERWISE 

25  LET  COM3IN=N 

26  LET  FACT =1.0 

27  FOR  Ir2  TO  N,  LET  FACT s FACT* I  ’’FOR  N  FACTORIAL 

23  LET  PDFsT**(N-1) 

29  LET  CDF=PDr*T 

30  LET  SIGSr  -1.0 

31  FOR  1=1  TO  N  DO 


A-2 


32  LET  TIrl 

33  IF  T  LE  TI 

3*4  GO  TO  LI 

35  OTHERWISE 

36  LET  TERM:SIGt)»COMBIN»(T-TI)»»(N-l) 

37  ADD  TERM  TO  PDF 

38  ADD  TERM»(T-TI)  TO  CDF 

39  LET  SIGN:  -SIGN 

40  LET  C0MBIN=C0MBIN»(N-I)/(I+1) 

41  LOOP  "OVER  I 

42  'Ll 'LET  PDFrPDF»N/FACT 

43  LET  CDFzCDF/FACT 

44  RETURN 

45  END  "NFOLD.U 

1  FUNCTION  SNORM(Z) 

2  " 

3  "ROUTINE  CALCULATES  THE  STANDARD  NORMAL  PROBABILITY  INTEGRAL. 

4  "REF:  APPROXIMATION  OBTAINED  FROM  AMS  55,  ABRAMOWITZ  AND  STEGUN. 

5  " 

6  IF  ABS.F(Z)  >  7.0 

7  GO  TO  L2 

8  OTHERWISE 

9  LET  P:0.5*SIGN.F(Z)»0.5»ERRFX(ABS.F(Z)/SQRT.F(2.0)) 

10  RETURN  WITH  P 

11  'L2'LET  P=0.5*SIGN.F(Z)»0.5 

12  RETURN  WITH  P 

13  END  "OF  SNORM 

1  FUNCTION  ERRFX(X) 

2  " 

3  "ROUTINE  CALCULATES  THE  ERROR  FU.NCTION.  THIS  FUNCTION  IS  CALLED  BY 

4  "SNORM(Z). 

5  "REFERENCE:  AMS  55,  'HANDBOOK  OF  MATHEMATICAL  FUNCTIONS',  NAT.  BUREAU 

6  "OF  STANDARDS,  NOV.  1970,  (P.  299). 

7  " 

8  LET  S:SIGN.F(X) 

9  LET  X:ABS.F(X) 

10  IF  X< 0.00000000001 

11  RETURN  WITH  0.0 

12  OTHERWISE 

13  IF  X>10.0 

14  RETURN  WITH  S 

15  OTHERWISE 

16  LET  T:1. 0/(1. 0*0. 327591l«X) 

17  LET  SUM:1.06140543»T 

18  LET  SUM: (SUM- 1.4531 5203 )»T 

19  LET  SUM: (SUM* 1.421 41 374 )»T 

20  LET  S'JM:(5UM-0.284496736)»T 

21  LET  SUM:(SUM*0.254829592)»T 

22  RETURN  WITH  S»(  1 . 0-SUM«EXP.F(-X''X) ) 

23  END  "OF  FUNCTION  ERRFX 
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ANNEX  B 


SIMSCRIPT  SOURCE  PROGRAM;  RUN.NFOLD.GU 

1  PREAMBLE  "RUN.NFOLD.GU 

2  NORMALLY  MODE  IS  REAL 

3  DEFINE  SNORM  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

4  DEFINE  ERRFX  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

5  DEFINE  ICOMBIN  AS  AN  INTEGER  FUNCTION  GIVEN  2  ARGUMENTS 

6  END  "PREAMBLE 


1  MAIN  "RUN.NFOLD.GU 

2  " 

3  "Driver  program  to  obtain  the  probability  density  and  cum  prob  of  the 

4  "N-fold  convolution  of  a  set  of  uniform  distributions  having  a  common 

5  "lower  domain  limit  (CL)  and  having  different  upper  domain  limits. 

6  "PDFs  &  CDFs  of  Normal  dist  having  same  avg  and  s.d.  are  also  printed. 

7  " 


8 

DEFINE  FLAGM,I.J,K,M,N,NCELLS,NREPS,SEED  AS  INTEGER 

VARIABLES 

9 

DEFINE  ANSWER  AS  A  TEXT  VARIABLE 

10 

DEFINE  NCV.HISTV  AS  INTEGER,  l-DIMENSIONAL  ARRAYS 

n 

DEFINE  AV,XV,CDFV  AS  REAL,  1-DIMENSIONAL  ARRAYS 

12 

DEFINE  SM  AS  A  REAL,  2-DIMENSIONAL  ARRAY 

13 

LET  LINES. Vr 9999 

14 

LET  RT12=SQRT.F(12.0) 

15 

LET  NCELLSrIO 

16 

RESERVE  CDFV(»)  AS  NCELLS 

17  'L0*SKIP  1  LINE 

18 

PRINT  7  LINES  THUS 

This  program  calculates  and  prints  the  p.d.f.  and  c.d.f  of  the 

N-fold  con- 

volution 

of  a  set  of  N  uniform  probability  distributions,  each  of  which  is 

defined  on  its  own,  possibly,  unique  interval— CL  to  upper  limit.  Inputs  are 

Integer  N 

,  upper  c.d.f.  value  to  terminate  calculation  (PMAX), 

common  lower 

argument 

value  (CL),  and  N  upper  limits  of  the  uniform  ranges. 

Max  value  of  i 

permitted 

by  the  program  Is  20.  Optionally,  a  Monte-Carlo  histogram  can  be 

obtained. 

26 

SKIP  2  LINES 

27 

PRINT  1  LINE  THUS 

INPUT 

THE  VALUE  OF  N. 

29 

READ  N 

30 

LET  N=MIN.F(N,20) 

31 

RESERVE  AV(»)  AS  N 

32 

RESERVE  NCV(»)  AS  N 

33 

LET  NCV(1)=N 

34 

LET  NCV(N)=1 

35 

FOR  K=2  TO  N-l,  LET  NCV(K)=ICOMBIN(N,K) 

36  " 

37  " 

RESERVE  MATRIX  OF  N-TUPLE  SUMS  OF  AV(»). 

38  " 

39 

RESERVE  SM(*,»)  AS  N  BY  • 

40 

FOR  1:1  TO  N,  RESERVE  SM(I,*)  AS  NCV(I) 

41 

PRINT  1  LINE  THUS 

INPUT  THE  (MAX)  VALUE  OF  THE  CONVOLUTION  CDF  TO  TERMINATE  CALCULATIONS. 

43  READ  PMAX 

44  PRINT  1  LINE  THUS 

INPUT  THE  C0rM)N  VALUE  OF  THE  ARGUMENT  LOWER  LIMIT  (OR  THRESHOLD). 

46  READ  CL 
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J7  LET  AVOsO.O 

l8  LET  VAH=0.0 

t9  LET  AC0N=1.0 

>0  FOR  1=1  TO  N  DO 

n  PRINT  1  LINE  WITH  I 

>2  THUS 

INPUT  THE  UPPER  LIMIT  OF  THE  ARGUMENT  RANGE  FOR  UNIFORM  VARIABLE  #  •». 


READ  AU 

55  IF  AU  LE  CL 

56  PRINT  1  LINE  WITH  AU,CL 

57  THUS 

INPUT  ERROR.  UPPER  ARG  LIM  .  IS  LESS  THAN  LOWER  . 

59  STOP 

60  OTHERWISE 

61  LET  AV(I)=AU-CL 

62  LET  ACON=ACON*AV(I) 

63  ADD  0.5*AV(I)  TO  AVG 

69  ADD  AV{I)»»2/12.0  TO  VAR 

65  LOOP  "OVER  (I)  UNIFORM  COMPONENTS 

66  LET  C0ND=1.0/SQRT.F(2.0»PI.C«VAB)  "FOB  NORMAL  DENSITY  COEF 

67  LET  STDVsSQRT.F(VAR) 

68  LET  LIM* AVG  ♦  3.0»STDV 

69  LET  LIM*MIN.F(LIM,  2.0»AVG) 

70  LET  DELT*LIM/20.0 

71  PRINT  1  LINE  THUS 

DO  YOU  WANT  A  MONTE-CARLO  SIMULATION?  (YES  OR  NO). 

73  READ  ANSWER 

79  IF  SUBSTR.F( ANSWER, 1,1)  *  "Y* 

75  LET  FLAGMsI 

76  PRINT  1  LINE  THUS 

INPUT  THE  INDEX  (1  TO  9)  OF  THE  RANDOM  NUMBER  SEED. 

78  READ  SEED 

79  PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  OF  REPLICATIONS  WANTED. 

81  READ  NREPS 

82  PRINT  1  LINE  WITH  NREPS 

83  THUS 

A  Monte-Carlo  slnulatlon  of  roplicationa  tiaa  begun. 

85  LET  NCELLS*10 

86  LET  DELX=2.0*0ELT 

87  RESERVE  XV(«)  AS  NCELLS 

88  RESERVE  HISTV(*)  AS  NCELLS 

89  FOR  K*1  TO  NCELLS,  LET  HI5TV(K)<0 

90  LET  AVGXtO.O 

91  LET  VARXzO.O 

92  FOR  K*1  TO  NCELLS,  LET  XV(K).N«CL*K*OELX 

93  " 

99  "SIMULATE  FOR  NREPS  REPLICATIONS. 


FOR  1:1  TO  NREPS  DO 
LET  SUM>0.0 
FOR  J*1  TO  N  DO 

ADD  UNIFORM. F(CL,CL*AV(J), SEED)  TO  SUM 
LOOP  "OVER  J 
ADO  SUM  TO  AVGX 
ADO  SUM**2  TO  VARX 
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104  *• DETERMINE  CELL  OF  HISTOGRAM  AND  ADD  TO  CELL  COUNT. 


105 

106 

107 

108 
109 

no 

111 

112  *K2* 

113 

114 


FOR  K=1  TO  NCELLS  DO 
IF  SUM  LE  XV(K) 

ADD  1  TO  HISTV(K) 
GO  TO  K2 
OTHERWISE 
LOOP  ’’OVER  K 

LOOP  "OVER  (I)  REPLICATIONS 

LET  AVGXrAVGX/NREPS 

LET  VARX=VARX/NREPS-AVGX«»2 


115  PRINT  1  LINE  THUS 

Monte-Carlo  simulation  has  been  completed. 

117  OTHERWISE 

118  LET  FLAGM=0 

119  ALWAYS 

120  " 

121  "FILL  RAGGED  TABLE  SM  WITH  N-TUPLE  SUMS  OF  THE  ELEMENTS  OF  AV. 


122  " 

123  CALL  STUPLES  (N,  AV(»),  SM(»,«)) 

124  SKIP  2  LINES 

125  PRINT  7  LINES  WITH  N 

126  THUS 


PROB  DISTRIBUTION  OF  A  ••-FOLD  CONVOLUTION  OF  A  SET  OF  UNIFORM  FUNCTIONS 


Indep  N- fold  Con vol ut ion  Normal  Prob  Distrib  *  Difference 

Variable  p.d.f.  c.d.f.  p.d.f.  c.d.f.  c.d.f. 


134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 


LET  RMS.DIFFsO.O 
LET  MAE.DIFFsO.O 
LET  KzO  "TO  COUNT  PAIRS 
FOR  Isl  TO  20  DO 
LET  TsI«DELT 

CALL  NFOLD.GU  (N,  ACON,  SH(^,*),  T)  YIELDING  PDF,  CDF 
IF  M0D.F(I,2)=0 
ADD  1  TO  K 
LET  CDFV(K)sCOF 
ALWAYS 

LET  ARG=(T-AVG)/STDV 

LET  NPDF=C0ND^EXP.F(-0.5*ARG^^2) 

LET  NCDF=SNORM(ARG) 

LET  DIfF=NCOF-CDF 

ADD  DIFF^^2  TO  RMS.DIFF 

LET  MAE.DIFFsMAX.F(MAE.DIFF,ABS.F(DIFF)) 

PRINT  1  LINE  WITH  T^N^CL,  PDF,  CDF,  NPDF,  NCDF,  DIFF 
THUS 


153 

154 

155 

156 

157 


'Ll 


».••••••  •.•••••• 

IF  CDF  GE  PMAX 
GO  TO  LI 
OTHERWISE 
LOOP  "OVER  I 
'PRINT  2  LINES  THUS 


160 


LET  RMS, DIFF=SgHT.F(0.05^RMS. DIFF) 
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161  PRINT  2  LINES  WITH  MAE.DIFF,RMS.DIFF 

162  THUS 

Max  abs  error  between  c.d.f.  and  the  Normal  c.d.f.  approx  »,•»»•»» 

RMS  difference  between  c.d.f,  and  the  Normal  c.d.f.  approx 

165  PRINT  3  LINES  WITH  AVG+N»CL,STDV 

166  THUS 

Mean  value  of  the  convolution  is  with  std  dev 

(1)  Mean  and  Std  Dev  of  each  of  the  Uniform  distributions: 

Component  Mean  Value  Std  Deviation 

170  FOR  1=1  TO  N  DO 

171  PRINT  1  LINE  WITH  I,  CL+0.5*AV(I) ,  AV(I)/RT12 

•  172  THUS 

*• 

174  LOOP  **0VER  (I)  UNIFORM  DISTRIBUTIONS 

175  SKIP  2  LINES 

176  IF  FLAGM  NE  1 

177  GO  TO  K3 

178  OTHERWISE 

179  PRINT  7  LINES  WITH  N.NREPS 

180  THUS 

MONTE-CARLO  SAMPLE  DISTRIBUTION  OF  THE  SUM  OF  ••  UNIFORM  RANDOM  VARIABLES 
NUMBER  OF  REPLICATIONS:  •••••• 


Indep 

Variable 

Histo  Sample  Sample  Dil*^  Versus 

Frequency  p.d.f.  c.d.f.  analytic  c.d.f. 

TOT 

LET 

XCUFsU.U 

189 

LET 

RMS.OIFFsO.O 

190 

LET 

MsO 

191 

FOR 

K.1  TO  NCELLS,  ADD  HISTV(K)  TO  M 

192 

FOR 

K«1  TO  NCELLS  DO 

193 

LET  XPDFxHISTV(K)/M 

194 

LET  XCDFsXCDF^XPDF 

195 

LET  DIFFsXCDF-COFV(K) 

196 

ADD  DIFF»»2  TO  RMS.DIFF 

197 

PRINT  1  LINE  NITH  XV(K),HISTV(K).XPOF,XCDF,DIFF 

198 

THUS 

200 

LOOP  "OVER  (K)  HISTO  CELLS 

201 

PRINT  2  LINES  THUS 

204  LET  RMS.OIFP=SQRT.F(RHS.OIFF/REAL.P(NCBLLS)) 

205  PRINT  1  LINE  WITH  RMS.DIFF 

206  THUS 

RMS  difference:  sample  c.d.f.  -  analytic  c.d.f.  e.eeeeee 

208  LET  SDX=SQRT.F(VARX) 

209  LET  SEX=SDX/SQRT.F(REAL.F(NREPS)) 

210  PRINT  3  LINES  WITH  AVGX,SDX»AVGX.1 .96*S£X,AVGX^1 .96«SEX 

211  THUS 

Sample  Average  Value  •••■•. eaeee  Sample  Standard  Deviation  eeeee^ei 
95  percent  confidence  interval  in  mean:  eeeee.eeeee^  eeeee,eeeee 

215  *K3*RELEASE  NCV(«) 

216  RELEASE  AV(«) 

217  RELEASE  SM(»,») 
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218  PRINT  1  LINE  THUS 

DO  YOU  HAVE  OTHER  PROBLEMS  OF  THIS  KIND?  (YES  OR  NO), 

220  READ  ANSWER 

221  IF  SUBSTR.FC ANSWER, 1,1)  =  «Y" 

222  GO  TO  LO 

223  OTHERWISE 

224  STOP 

225  END  ’'MAIN 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 


ROUTINE  NFOLD.GU  GIVEN  N,  AGON,  SM,  T  YIELDING  PL/,  CDF 

» I 


Routine  calculates  the  prob  density  function  (PDF)  and  the  cum¬ 
ulative  dlst  function  (CDF)  of  the  N-fold  convolution  of  a  set  of 
N  unique  uniform  distributions.  The  range  of  the  1  th  distribution 
"Is  (0,  a(l)).  The  product,  over  N,  of  the  a(l)  Is  the  argument  ACON 
All  n-tuple  sums  of  elements  a(l)  are  entered  In  the  ragged  table  SM 
where  the  k  th  row  and  j  th  column  element  Is  the  J  th  k- tuple  sum. 
E.g.,  the  first  row  of  SM  contains  a(J),  Real- valued  argument  Is  T. 
Num  of  combinations  of  N  objects  taken  K  at  a  time  Is  DIM.F(SM(K,*)) 
With  the  following  notation  for  the  CDF,  with  argument  tt  F(n,t), 
with  the  j  th  k-tuple  sum  for  the  n  th  convolution  denoted  by 


S  (n) 
kj 


and  with  the  unit  step  function  at  x  denoted  by  u(t-x). 


F(n,t)  s  (1/AC0N/n!)(t  Sum  over  kt1  to  n  and  js1  to  C(n,k): 


(-1)  u(t  -  S  (n))(t  - 
jk 


n 

S  (n))  ), 
jk 


where  C(n,k)  Is  the  #  combinations  of  n  things  taken  k  at  a  time. 


DEFINE  I,J,K,N  AS  INTEGER  VARIABLES 
DEFINE  SM  AS  A  REAL,  2-DIMENSIONAL  ARRAY 
IF  T  LE  0.0 
LET  PDFsO.O 
LET  CDFsO.O 
RETURN 
OTHERWISE 
IF  T  GE  SM(N,1) 

LET  PDFrO.O 
LET  CDFrI.O 
RETURN 
OTHERWISE 
LET  FACTS  1,0 

FOR  1=2  TO  N,  LET  FACT=FACT«I  "FOR  N  FACTORIAL 
LET  PDF=T»»(N-1) 

LET  C0FsP0F»T 
LET  SIGNS  1.0 
FOR  K=1  TO  N  DO 
LET  SIGNS  -SIGN 
FOR  J=1  TO  0IM.F(SM(K,»))  DO 
IF  T  >  SM(K,J) 
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U9 

LET  TARGsT-SM(K,J) 

50 

LET  INCR=SION*TARG*«N 

51 

ADD  INCR/TARG  TO  PDF 

52 

ADD  INCR  TO  CDF 

53 

ALWAYS 

54 

LOOP  "OVER  (J)  COLUMNS 

55 

LOOP  "OVER  (K)  ROWS 

56 

LET  PDFsPDF»N/FACT/ACON 

57 

LET  CDFrCDF/FACT/ACON 

58 

RETURN 

59 

END  "NFOLD.GU 

1 

5 

FUNCTION  ICONBIN  (N,  K) 

1  1 

C 

3 

li 

*  *  INTEGER-VALUED  #  OP  COMBINATIONS  OP  N  OBJECTS 

1 1 

TAKEN  K  AT  A  TIME. 

5 

DEPINE  C,I«K,N  AS  INTEGER  VARIABLES 

6 

IP  X  «  0 

7 

RETURN  WITH  1 

8 

OTHERWISE 

9 

LET  C*1 

10 

POR  Is1  TO  K  DO 

11 

LET  C*C»(N.l4l)/I 

12 

LOOP  ‘‘OVER  I 

13 

RETURN  WITH  C 

14 

END  ••PUNCTION  ICOMBIN 

1 

ROUTINE  STUPLES  (N.  AV,  SH) 

t  V 

C 

3 

*' Routine  fills  the  slesents  of  a  ragged  tablt» 

SH.  having  N  rows* 

u 

"The  altaant  of  this  tabla  consists  of  the  J  th  k-tupla  sua  of 

5 

"•laments  of  the  vector  AV*  Routine  is  called 

V  1 

by  NFOLD.GU. 

0 

7 

DEFINE  1, 11, 12.13, IR. 15.16,17. X8.I9iI10.in, 

112.113,11*1.115,116 

8 

I17.I18.I19.J,K,N  AS  INTEGER  VARIABLES 

9 

DEPINE  JV  AS  AN  INTEGER.  ^DIMENSIONAL  ARRAY 

10 

DEPINE  AV  AS  A  REAL.  I-DIMENSIONAL  ARRAY 

11 

DEFINE  SM  AS  A  REAL.  2-DIHENSIONAL  ARRAY 

12 

IP  N  >  20 

13 

PRINT  1  LINE  WITH  N 

14 

THUS 

ilUHBER  OF  VARIABLES  («  ••)  EXCEEDS  THE  CAPACITY  OF  20  IN  ROUTINE  STUPLES. 

16 

STOP 

17 

OTHERWISE 

18 

RESERVE  JV(»)  AS  N  "LOCAaY 

19 

LET  SM(N.1)«0.0 

20 

POR  Jf1  TO  N.  ADO  AV(J}  TO  SHCN.I) 

21 

FOR  Iltl  TO  N  DO 

22 

LET  SItAVdl) 

23 

LET  SMd.ID.SI 

24 

IP  N  <  3 

25 

CO  TO  LI 

26 

OTHERWISE  "gen  2  tuples 

27 

POR  I2sIU1  TO  N  DO 

28 

ADO  1  TO  JV(2) 

29 

LET  S2*SUAV(I2) 

30 

LET  SMC2.JV(2))sS2 

,1 


31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 
87 


IF  N  <  4 
GO  TO  L2 

OTHERWISE  ”gen  3  tuples 
FOR  13=12+1  TO  N  DO 
ADD  1  TO  JV(3) 

LET  S3=S2+AV(I3) 

LET  SM(3,JV(3))=S3 
IF  N  <  5 
GO  TO  L3 

OTHERWISE  "gen  4  tuples 
FOR  I4=I3-*»1  TO  N  DO 
ADD  1  TO  JV(4) 

LET  S4=S3+AV(I4) 

LET  SM(4,JV(4))=S4 
IF  N  <  6 
GO  TO  L4 

OTHERWISE  "gen  5  tuples 
FOR  15=14+1  TO  N  DO 
ADD  1  TO  JV(5) 

LET  S5=S4+AV(I5) 

LET  SM(5,JV(5))=S5 
IF  N  <  7 
GO  TO  L5 

OTHERWISE  "gen  6  tuples 
FOR  I6=I5-^1  TO  N  DO 
ADD  1  TO  JV(6) 

LET  S6=S5^»AV(I6) 

LET  SM(6,JV(6))=S6 
IF  N  <  8 
GO  TO  L6 

OTHERWISE  "gen  7  tuples 
FOR  17=16+1  TO  N  DO 
ADD  1  TO  JV(7) 

LET  S7=S6+AV(I7) 

LET  SM(7,JV(7))=S7 
IF  N  <  9 
GO  TO  L7 

OTHERWISE  "gen  8  tuples 
FOR  I8=I7*1  TO  N  DO 
ADD  1  TO  JV(8) 

LET  S8=S7+AV(I8) 

LET  SM(8,JV(8))=S8 
IF  N  <  10 
GO  TO  L8 

OTHERWISE  "gen  9  tuples 
FOR  19=18+1  TO  N  DO 
ADD  1  TO  JVC 9) 

LET  S9=S8+AV(I9) 

LET  SM(9,JV(9))=S9 
IF  N  <  11 
GO  TO  L9 

OTHERWISE  "gen  10  tuples 
FOR  110=19+1  TO  N  DO 
ADD  1  TO  JVC  10) 

LET  S10=S9+AVC110) 

LET  SHC10,JVC10))=S10 
IF  N  <  12 
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88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 
103 
10U 

105 

106 

107 

108 
109 
no 
in 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 


GO  TO  L10 

OTHERWISE  '*gen  11  tuples 
FOR  111=110+1  TO  N  DO 
ADD  1  TO  JVC  11) 

LET  S11=S10+AV(I11) 

LET  SM(11,JV(11))=S11 
IF  N  <  13 
GO  TO  L11 

OTHERWISE  ‘'gen  12  tuples 
FOR  112=111+1  TO  N  DO 
ADD  1  TO  JVC  12) 

LET  S12=S11+AVCI12) 

LET  SMC12,JVC12))=S12 
IF  N  <  14 
GO  TO  L12 

OTHERWISE  '»gen  13  tuples 
FOR  113=112+1  TO  N  DO 
ADD  1  TO  JVC  13) 

LET  S13=S12+AVCI13) 

LET  SMC13,JVC13))=S13 
IF  N  <  15 
GO  TO  LI 3 

OTHERWISE  “gen  14  tuples 
FOR  114=113+1  TO  N  DO 
ADD  1  TO  JVC  14) 

LET  S14=S13+AVCI14) 

LET  SMC14,JVC14))=S14 
IF  N  <  16 
GO  TO  LI 4 

OTHERWISE  “gen  15  tuples 
FOR  115=114+1  TO  N  DO 
ADD  1  TO  JVC  15) 

LET  S15=S14+AVCI15) 

LET  SMC15,JVC15))=S15 
IF  N  <  17 
GO  TO  L15 

OTHERWISE  “gen  16  tuples 
FOR  116=115+1  TO  N  DO 
ADD  1  TO  JVC  16) 

LET  S16=S15+AVCI16) 

LET  SMC16,JVC16))=S16 
IF  N  <  18 
GO  TO  L16 

OTHERWISE  “gen  17  tuples 
FOR  117=116+1  TO  N  DO 
ADD  1  TO  JVC  17) 

LET  S17*^S16+AVCI17) 

LET  SMC17,JVC17))=S17 
IF  N  <  19 
GO  TO  L17 

OTHERWISE  “gen  18  tuples 
FOR  118=117+1  TO  N  DO 
ADD  1  TO  JVC  18) 

LET  S18=S17+AVCI18) 

LET  SMC18,JVC18))=S18 
IF  N  <  20 
GO  TO  L18 
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i)) ' 


145 

OTHERWISE  "gen  19  tuples 

146 

FOR  119=118+1  TO  N  DO 

147 

ADD  1  TO  JVC  19) 

148 

LET  S19=S18+AV(I19) 

149 

LET  SM(19,JV(19))=S19 

150 

•L19 

LOOP  "OVER  119 

151 

•L18 

LOOP  "OVER  118 

152 

’L17 

LOOP  "OVER  117 

153 

'L16 

LOOP  "OVER  116 

154 

•L15 

LOOP  "OVER  115 

155 

•L14 

LOOP  "OVER  114 

156 

•L13 

LOOP  "OVER  113 

157 

•L12 

LOOP  "OVER  112 

158 

'L11 

LOOP  "OVER  111 

159 

•L10 

LOOP  "OVER  110 

160 

•L9* 

LOOP  "OVER  19 

161 

'L8* 

LOOP  "OVER  18 

162 

•L7’ 

LOOP  "OVER  17 

163 

•L6‘ 

LOOP  "OVER  16 

164 

’L5’ 

LOOP  "OVER  15 

165 

•L4' 

LOOP  "OVER  14 

166 

'L3' 

LOOP  "OVER  13 

167 

'L2' 

LOOP  "OVER  12 

168 

’LTI 

LOOP  "OVER  11 

169 

RELEASE  JV(*) 

170 

RETURN 

171 

END 

'  'STUPLES 

1  FUNCTION  SNORM(Z) 

2  »i 

3  "ROUTINE  CALCULATES  THE  STANDARD  NORMAL  PROBABILITY  INTEGRAL. 

4  "REF:  APPROXIMATION  OBTAINED  FROM  AMS  55,  ABRAMOWITZ  AND  STEGUN. 

5  " 

6  IF  ABS.F(Z)  >  7.0 

7  GO  TO  L2 

8  OTHERWISE 

9  LET  P=0.5+SIGN.F(Z)»0.5»ERRFX(ABS.F(Z)/SQRT.F(2.0)) 

10  RETURN  WITH  P 

11  'L2'LET  P=0.5+SIGN.F(Z)«0.5 

12  RETURN  WITH  P 

13  END  "  OF  SNORM 


ANNEX  C 


SIMSCRIPT  SOURCE  PROGRAM:  RUN.NFOLD.E 

1  PREAMBLE  ’’RUN.NFOLD.E 

2  NORMALLY  MODE  IS  REAL 

3  DEFINE  SNORM  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

4  DEFINE  ERRFX  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

5  END  "PREAMBLE 

1  MAIN  "RUN.NFOLD.E 

2  " 

3  "Driver  program  to  obtain  the  probability  density  and  cam  prob  of  the 

4  "N-fold  convolution  of  a  set  of  N  exponential  dist's.  Two  options 

5  "are  available:  (a)  all  exponential  dist’s  in  the  set  to  be  convolved 

6  ’’are  identical,  (b)  all  of  the  exponential  dist's  are  unique.  PDFs 

7  ' ' and  CDFs  of  a  Normal  dist  with  the  same  mean  and  variance  is  printed 

8  ’’for  comparison  with  the  convolution. 

9  " 

10  DEFINE  I,J,K,L,M,N,NCELLS  AS  INTEGER  VARIABLES 

11  DEFINE  ANSWER  AS  A  TEXT  VARIABLE 

12  DEFINE  INDEX  AS  AN  INTEGER,  1-DIMENSIONAL  ARRAY 

13  DEFINE  AV,LAV  AS  REAL,  1-DIMSNSIONAL  ARRAYS 

14  DEFINE  MAT  AS  A  REAL,  2.DIMENSI0NAL  ARRAY 

15  LET  MCELLS=20 

16  ’LO'SKIP  1  LINE 

17  PRINT  7  LINES  THUS 

This  progr’am  calculates  and  prints  the  p.d.f.  and  c.d.f  of  the  N-fold 
convolution  of  a  set  of  N  exponential  distributions.  Two  progam  options  exist 
for  these  distributions:  (a)  all  distributions  have  the  same  mean  value,  and 
(b)  all  distributions  have  unique  (or  different)  means.  User  inputs  are  tne 
number  (N)  of  distributions  to  be  convolved,  the  upper  c.d.f,  limit  to  term¬ 
inate  calculations ,  and  the  mean  values  of  these  exponential  distributions, 

25  PRINT  1  LINE  THUS 

INPUT  THE  VALUE  0?  N, 

27  READ  N 

28  PRINT  1  LINE  THUS 

INPUT  THE  (MAX)  VALUE  OF  THE  CDF  TO  TERMINATE  CALCULATIONS. 

30  READ  PMAX 

31  RESERVE  AV(<*),LAV(‘»)  AS  N 

32  RESERVE  INDEX(»)  AS  N-1 

33  RESERVE  MAT(»,*)  AS  N  BY  N 

34  PRINT  1  LINE  THUS 

DO  ALL  EXPONENTIAL  DISTRIBUTIONS  HAVE  THE  SAME  PARAMETER?  (Y  OR  N). 

36  READ  ANSWER 

37  IF  SUBSTR.FC ANSWER, 1,1)  =  "Y" 

38  PRINT  1  LINE  THUS 
INPUT  THE  CO:WOM  MEAN  VALUE, 


40 

READ  AVG 

41 

FOR 

1=1  TO  N,  LET  LAV(I)=1.0/AVG 

42 

LET 

VAR  =  :]*»AVG’*»2 

4i 

LET 

AVG=N«AVG 

4-4 

OTHERWISE 

45 

LET 

AVGrO.O 

46 

f  tTT 

VARrO.O 

47 

LET 

LAMBDA =1.0 
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43  IF  N  >  15 

49  PRIWT  1  LINE  THUS 

Mx\X  0  OF  CONVOLUTIONS  IS  LIMITED  TO  15  FOR  THIS  OPTION. 

51  RELEASE  AV(»)  ,LAV(»)  ,INDSX(’») 

52  RELEASE  MAT(^,») 

53  LET  N=15 

54  RESERVE  AV(») ,LAV(») ^INDEXC*)  AS  N 

55  RESERVE  MAT(*,»)  AS  N  BY  N 

56  ALWAYS 

57  FOR  1=1  TO  N  DO 

53  PRINT  1  LINE  WITH  I 

59  THUS 

INPUT  THE  MEAN  VALUE  OF  THE  TH  EXPONENTIAL  DISTRIBUTION. 


CHOOSE  THE  15. 


61 

62 

63 

64 

65 

66 
67 
63 

69 

70 

71 

72 

73 

74 

75 

76 

77 
73 

79 

80 
81 
82 

83 

84 

85 
80 
87 
S8 

89 

90 

91 

92 

93 

94 

95 

96 

97 
Jib 
9') 

100 


READ  THETA 
ADD  THETA  TO  AVG 
ADD  THETA»»2  TO  VAR 
LET  LAV(I)=1.0/THETA 
LET  LAMBDA=LAMBDA»LAV(I) 

LOOP  ’’OVER  (I)  EXPONENTIAL  DISTRIBUTIONS 

1 1 

"OBTAIN  THE  COEFFICIENTS  AV(»)  FOR  THE  CONVOLUTION  DENSITY. 

t  » 

IF  N  =  2 

LET  AV(1)=LAMBDA/(LAV(2)-LAV(1)) 

LET  AV(2)=  -AV(1) 

00  TO  L2 
OTHERNISE 

FOR  J=1  TO  N,  LET  MAT(1,J)=1.0 
FOR  1=2  TO  N-1,  FOR  J=1  TO  N,  LET  !<AT(I, J)=0.0 
FOR  J=1  TO  N  DO 
LET  PNJsl.O 
FOR  Ksl  TO  N  DO 
IF  K  NB  J 

LET  PNJ=PNJ.LAV(K) 

ALWAYS 

LOOP  "OVER  K 
LET  MAT(N,J)=PNJ 
LOOP  "OVER  (J)  COLUMNS 

CALL  HTUPLES  (N,  LAVC*),  MAK*,"*))  "FILL  EL'MTS  OF  MAT(»,») 

I  I 

"OBTAIN  THE  INVERSE  OF  MAT («,»). 

•  I 

CALL  MAT. INVERSE  (N,  MAT(»,»)) 

FOR  1=1  TO  N,  LET  AV(I)=LAMBDA»MAT(I ,N) 

ALWAYS 

’ L2 ’ LET  COND=  1 . 0/SORT . F( 2. 0<»PI .C^VAR ) 

LET  S7DV=SQRT.F(VAR) 

LET  LIMrAVG  ♦  4.0»STDV 
LET  DELTrLIM/NCELLS 
LET  LINES. V=9999 
SKIP  2  LINES 
PRINT  7  LINES  WITH  N 
THUS 


(1) 


PRJH  DISTRIBUTION  0?  THE  CONVOLUTION  OF  A  SET  OF  »»  EXPONENTIAL  DISTRIBUTIONS 
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IndTap”  N-f^old  Convorut ion  Norinai  Prob  Distrib  Difference 

Variable  p.d.f.  c.d.f.  p.d.f.  c.d.f.  c.d.f. 

f08  LEf’M!\E="oro 

109  LET  RMSrO.O 

110  FOR  1=1  TO  NCELLS  DO 

111  LET  T=I*^DELT 

112  CALL  NFOLD.E  (N,  AV(«),  LAV(»),  T)  YIELDING  PDF, CDF 

113  LET  ARG=(T-AVG)/STDV 

114  LET  NPDF=C0ND^*EXP.F(-0.5»ARG**2) 

115  LET  NCDF=SNORM(ARG) 

116  LET  DIFF=CDF-NCDF 

117  LET  MAE=MAX.F(MAE,ABS.F(DIFF)) 

118  ADD  DIFF»»2  TO  RMS 

119  PRINT  1  LINE  WITH  T,  PDF,  CDF,  NPDF,  MCDF,  DIFF 

120  THUS 

122  IF  CDF  GE  PMAX 

123  GO  TO  LI 

124  OTHERWISE 

125  LOOP  "OVER  I 

126  ’Ll ’PRINT  2  LINES  THUS 


129  LET  RMS=SQRT.F(RMS/REAL.F(NCELLS)) 

130  PRINT  2  LINES  WITH  tlAE,RMS 

131  THUS 

Max  abs  error  In  a  Normal  approximation  of  sum  c.d.f.  • 
RMS  error  of  a  Normal  approximation  of  the  sum  c.d.f.  ^  * 

134  PRINT  3  LINES  WITH  AVG,  STDV 

135  THUS 

(1)  Mean  Value  of  the  Convolution  .  Std  Dev 

Mean  Values  of  each  of  the  exponential  distributions 


139  FOR  1=1  TO  N  DO 

140  PRINT  1  LINE  WITH  I,  1.0/LAV(I) 

141  THUS 

Number  **  Mean  . 

143  LOOP  "OVER  I 

144  SKIP  2  LINES 

145  RELEASE  LAV(») 

146  RELEASE  AV(») 

147  RELEASE  MAT(»,») 

148  PRINT  1  LINE  THUS 

DO  YOU  HAVE  SIMILAR  PROBLEtiS  TO  SOLVE?  (YES  OR  NO). 

150  READ  ANSWER 

151  IF  SUBSTR.FC ANSWER, 1,1)  =  "Y” 

152  GO  TO  LO 

153  OTHERWISE 

154  STOP 

155  FND  "MAIN 
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1 

2 

3 

n 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 


ROUTINE  .\TOLD.E  GIVEM  N,  AV,  LAV,  T  YIELDING  PDF,  CDF 


20 

21 

2^ 

23 

24 
23 
26 
27 


29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 
49 

49 

50 


52 

53 


Routine  calculates  the  probability  density  function  (PDF)  and  the  can 
distribution  function  (CDF)  of  the  N-fold  convolution  of  a  set  of 
exponential  dist's.  Two  options  are  provided:  (a)  each  of  the  N 
dist's  nas  the  sane  mean,  and  (b)  each  dist  has  a  unique  mean. 

If  the  mean  values  of  the  exponential  dist's  are  all  unique,  the 
convolution  can  be  expressed  as  a  weighted  san  of  the  exponential 
dist's.  These  weignts  are  passed  to  the  routine  in  the  vector  AV(*). 
Rate  parameters  of  the  N  expon  dists  are  given  in  the  vector  LAV(*). 
The  real- valued  argument  of  the  convolution  is  T. 


If  each  of  the  dist's  has  the  same  rate  parameter,  a,  (option  (a)), 
tne  p.d.f,  and  c.d.f.  of  the  n-foid  convolution  are  given, 
respectively,  by  these  Erlang  functions: 


END 


n- 1 

f(n,t)  =  a  (at)  exp(-at )/ (n- 1 ) ! 


t  >  0. 


F(n,t)  =  1  -  exp(-at)  SaTi(i=0  to  n-1):  z  /  i! 


DEFINE  I,N  AS  INTEGER  VARIABLES 

DEFINE  AV,LAV  AS  REAL,  1-DIMENSIONAL  ARRAYS 

IF  T  LE  0.0 

LET  PDFrO.O 
LET  CDF=0.0 
RETURN 
OTHERWISE 

IF  LAV(1)=LAV(2)  '’ALL  RATE  PARMS  ASSUMED  EQUAL 
LET  Z=LAV(l)n 
LET  EXP2=EXP.F(-Z) 

LET  FACT =1,0 
LET  ZIzl.O 
LET  SUMzI.O 
FOR  1=1  TO  N-1  DO 
LET  FACT=FACT»I 
LET  ZI=Z»ZI 
ADD  ZI/FACT  TO  SUM 
LOOP  "OVER  I 

LET  PDF=LAV(1)*ZI/FACT»EXPZ 
LET  CDF=1,0-EXPZ»SUM 
RETURN 
OTHERWISE 
'.ET  PDFrO.O 
LET  CDFrO.O 
FOR  1=1  TO  N  DO 

LET  EXP:=EXP.F(-LAV(I)*T) 

ADD  AV(I)»EXPZ  to  PDF 
ADD  AV(I)/LAV(I)»(1.0-EXPZ)  TO  CDF 
LOOP  "OVER  I 
RETURN 
"NFOLD.E 
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1  ROUTINE  FOR  MAT. INVERSE  (N,  AM) 

2  ’  ’ 

3  ’’ROUTINE  TO  OBTAIN  THE  INVERSE  0?  THE  N  BY  N  MATRIX  AM  VIA  THE 

4  ’’COMPACT  FORM  OF  THE  GAUSS-JORDAN  METHOD.  INV  IS  RETURNED  IN  AM. 

5  ” 

6  DEFINE  I,  J,  K,  M  AS  INTEGER  VARIABLES 

7  DEFINE  AM  AS  A  REAL,  2-DIMEMSIONAL  ARRAY 

8  FOR  1=1  TO  N  DO 

9  LET  P=AM(I,I) 

10  IF  P=0.0 

11  PRINT  2  LINES  WITH  I  THUS 

ERROR  IN  ROUTINE  MAT. INVERSE.  THE  TH  DIAGONAL  ELEMENT  IS  ZERO. 

THE  MATRIX  CANNOT  BE  INVERTED. 


14 

STOP 

15 

OTHERWISE 

16 

LET  AM(I,I):1.0 

17 

FOR  J=1  TO  >)  DO 

18 

LET  AM(I,J)=AM(I,J)/P 

19 

LOOP  "OVER  J 

20 

FOR  J=1  TO  N  DO  ’’THE  SECOND  J-LOOP 

21 

IF  J  =  I 

22 

GO  TO  EOJ  ’ ’ END  0? 

J-LOOP 

23 

OTHERWISE 

24 

LET  P=AM(J,I) 

25 

LET  AM(J,I)=0.0 

26 

FOR  K=1  TO  N  DO 

27 

SUBTRACT  P»AM(I,K) 

FROM  AIi(J 

28 

LOOP  ’’OVER  K 

29 

’EOJ’  LOOP  ’’OVER  J 

30 

LOOP  ’’OVER  I 

31 

RETURN 

32 

END  ’’ROUTINE  MAT. INVERSE 

1 

ROUTINE  NTUPLES  (N,  LAV,  MAT) 

2  ” 

3  ’’Routine  fills  N-2  row  elements  in  the  MAT(*,**)  which  arc  contrlbuiOvj 

4  ’’by  n-tuples  associated  with  variable  L^V(**).  Routine  is  called  by 

5  ’’RUN.NFOLD.E. 

6  DEFINE  I, II, 12, 13, 14, 15, 16, 17, 18, 19, no, 111, 112, J,K,N  AS  INTEJEI: 
VARIABLES 

7  DEFINE  INDEX  AS  AN  INTEGER,  1-DIMENSIONAL  ARRAY 

3  DEFINE  LAV  AS  A  REAL,  1-DIMENSIONAL  ARRAY 

9  DEFINE  MAT  AS  A  REAL,  2-DIMSNSIONAL  ARRA^ 

10  IF  M  >  15 

11  PRINT  1  LINE  vJITH  N 

12  THUS 

INPUT  ERROR  TO  ROUTINE  NTUPLES.  NUMBER  OF  CON VOL’JT IONS,  »»,  IS  EXCESSIVE. 

1 4  STOP 

1 b  OTHERWISE 

16  RESERVE  INDEX{»)  AS  N-1  ’’LOCALLY 

17  FOR  J=1  TO  N  DO 

18  LET  K=0 

19  FOR  1=1  TO  N  DO 

20  IF  I  NE  J 

21  ADD  1  TO  X 

?.?  LET  INDEX(K)  =  I 

23  ALWAYS 


34 

35 
26 


28 

29 

30 

31 

32 

33 

34 

35 

36 

37 
33 

39 

40 

41 

42 

43 

44 

45 
40 
47 
4B 

49 

50 

51 

52 

53 

54 

55 


56 


57 

53 

59 

60 
61 
62 
63 
6*1 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 


LOOP  ’’OVEIl  (I)  PARAMETERS 
IF  N  <  3 

RELEASE  INDEX (»») 

RETURN 
OTHERWISE 
FOR  11=1  TO  N-1  DO 

LET  LA1=LAV(INDSX(I1)) 

ADD  LAI  TO  MAT(2,J) 

IF  N  <  4 
GO  TO  LI 

OTHERWISE  ”gen  2  tuples 
FOR  I2  =  IU1  TO  N-1  DO 

LET  LA2=LA1*LAV(INDEX(I2)) 

ADD  LA2  TO  MAT(3,J) 

IF  N  <  5 

30  TO  L2 

OTHERWISE  ’’sen  3  tuples 
FOR  13=12+1  TO  N-1  DO 

LET  LA3=LA2»LAV(INDSX(I3)) 

ADD  LAS  TO  MAT(4,J) 

IF  N  <  6 
GO  TO  L3 

OTHERWISE  "gen  4  tuples 
FOR  14=13+1  TO  N-1  DO 

LET  LA4=LA3**LAV(INDEX(I4)) 

ADD  LA4  TO  MAT(5,J) 

IF  N  <  7 
GO  TO  L4 

OTHERWISE  "gen  5  tuples 
FOR  15=14+1  TO  N-1  DO 

LET  LA5=U4»LAV(INDEX(I5)) 

ADD  LAS  TO  MAT(6,J) 

IF  N  <  8 
GO  TO  L5 

OTHERWISE  "gen  6  tuples 
FOR  16=15+1  TO  N-1  DO 

LET  LA6=LA5«LAV( INDEX (16)) 

ADD  LA6  TO  MAT(7,J) 

IF  N  <  9 
GO  TO  L6 

OTHERWISE  "gen  7  tuples 
FOR  17=16+1  TO  N-1  DO 

LET  LA7=LA6«LAV( INDEX (17)) 

ADD  LA7  TO  MAT(8,J) 

IF  N  <  10 
GO  TO  L7 

OTHERWISE  "gen  8  tuples 
FOR  18=17+1  TO  N-1  DO 

LET  LA8=LA7«LAV(INDEX(I8)) 

ADD  LA8  TO  MAT(9,J) 

IF  N  <  11 
GO  TO  L8 

OTHERWISE  "gon  9  tuples 
FOR  19=18+1  TO  N-1  DO 

LET  LA9=LA8*LAV(INDEX(I9)) 
ADD  LA9  TO  MAT(IO.J) 
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IF  N  <  12 
GO  TO  L9 

OTHERWISE  **gen  10  tuples 
FOR  110=19+1  TO  N-1  DO 

LET  LA 10=LA9*LAV( INDEX 
ADD  LA10  TO  MAT(11,J) 
IF  N  <  13 
GO  TO  L10 

OTHERWISE  "gen  11  tup 
FOR  111=110^.1  TO  N-1  D 
LET  LA11=LA10»LAV(INDE 
ADD  LA11  TO  MAT(12,J) 
IF  N  <  14 
GO  TO  L11 

OTHERWISE  "gen  12  tup 
FOR  I12=I1U1  TO  N-1  D 
LET  LA12=LA11^LAV(I 
ADD  LA12  TO  MAT(13, 

•LI 2’  LOOP  "OVER  112 

•Lir  LOOP  ••OVER  111 

•LIO^  LOOP  ••OVER  110 

•L9'  LOOP  ••OVER  19 

•L8^  LOOP  ••OVER  18 

•L7*  LOOP  ••OVER  17 

•L6^  LOOP  '•OVER  16 

•L5'  LOOP  ••OVER  15 

•L4^  LOOP  ••OVER  14 

•L3'  LOOP  ••OVER  13 

•L2»  LOOP  ••OVER  12 

•LI 'LOOP  "OVER  II 

LOOP  "OVER  (J)  COLU:^S 
RELEASE  INDEX (*) 

END  "ROUTINE  NTUPLES 
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SINSCfilPT  SOURCE  PROGRAM:  LP.INV 

1  PREAMBLE  LP.INV 

2  NORMALLY  NODE  IS  REAL 

3  DEFINE  A1,A2  AS  REAL  VARIABLES 

4  DEFINE  SNORM  AS  A  REf>L  FUNCTION  GIVEN  1  ARGUMENT 

5  DEFINE  ERRFX  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

6  DEFINE  LTRNS.FUN  AS  A  REAL  FUNCTION  GIVEN  3  ARGUMENTS 

7  END  ''PREAMBLE 

1  MAIN  "LP.INV 

2  " 

3  "Obtains  the  Inverse  Laplace  transform  at  a  set  of  discrete  points  via 

4  "Bellman's  method.  This  method  approximates  an  Integral  by  gausslan 

5  "quadrature.  The  quadrature  formula  leads  to  a  matrix  eq'n  whose 

6  "sol'n  defines  the  Inv  In  terms  of  the  transform  evaluated  at  M  points. 

7  "An  example  of  this  method  Is  provided  on  pp  14  thru  18  of  "Numerical 

8  "Methods  In  Renewal  Theory,"  (AD  828276),  Feb  1968. 

9  " 

10  DEFINE  FLAGE,FLAGN,FLAGMX,I,J,K,L,M,N,NREPS,SEED  AS  INTEGER 
VARIABLES 

11  DEFINE  ANSWER, FILNAM,TITLE,DIST.NAME  AS  TEXT  VARIABLES 

12  DEFINE  IPVT,HISTV  AS  INTEGER,  1-DIMENSIONAL  ARRAYS 

13  DEFINE  DET,TV  AS  REAL,  1-DIMENSIONAL  ARRAYS 

14  RESERVE  DET(*)  AS  2 

15  DEFINE  DSTARV,DV,GSTARV,GV,UV,WV,XV,CDFV  AS  REAL,  1-DIMENSIONAL 
ARRAYS 

16  DEFINE  AM  AS  A  REAL,  2-DIMENSIONAL  ARRAY 

17  LET  Ms  16  "TERMS  IN  THE  GAUSSIAN  QUADRATURE 

18  RESERVE  CDPV(»)  AS  M 

19  LET  FILNAH  s  "GAUSS. Q1 6. DATA" 

20  "  LET  DIST.NAME  s  "Gaaiu(3)" 

21  LET  DIST.NAME  s  "Expon  Mix" 

22  LET  FLAGMXsl 

23  LET  Ks1 

24  RESERVE  TV(*) ,HISTV(»)  AS  M 

25  LET  CON.AVGsK  "CONSTANT  RELATING  AVG  TO  RATE  FARM 

26  RESERVE  WV(»),  XV(*),  IPVT(»)  AS  M 

27  RESERVE  DSTARV(»),DV(»),GSTARV(«),GV(«)  AS  M 

28  RESERVE  AH(*,«)  AS  M  BY  H 

29  •* 

30  "READ  THE  QUADRATURE  POINTS  AND  WEIGHTS  FROM  THE  FILE:  FILNAH. 

31  •• 

32  LET  EOF. Vs  1 

33  LET  LINES. Vs9999 

34  OPEN  UNIT  4  FOR  INPUT, 

35  OLD, 

36  FILE  NAME  IS  FILNAH 

37  RECORD  SIZE  IS  120 

38  USE  UNIT  4  FOR  INPUT 

39  READ  TITLE  USING  UNIT  4 

40  PRINT  2  LINES  WITH  FILNAH, TITU 

41  THUS 

Data  file  eeeeeeeeetette  13  ^33^  for  eeeeeeeeeeeteeeeeeeeeeeeeeeeeeeeeeeeeee 
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««  FOB  1:1  TO  M  DO 

US  READ  XV(I),WV(I)  USING  UNIT  4 

46  LOOP  "OVER  (I)  QUADRATURE  POINTS 

47  CLOSE  UNIT  4 

48  USE  UNIT  5  FOR  INPUT 

49  FOR  1=1  TO  M,  LET  TV(I)=L0O.E.F(2.0/{XV{M-I+1)*1.0)) 

50  PRINT  7  LINES  WITH  DIST.NAME 

51  THUS 

This  program  calculates  and  prints  the  o.d.f .  of  an  N-fold  convolution  of 
a  set  of  N  distributions  with  different  mean  values.  This  o.d.f. 

Is  obtained  by  numerically  Inverting  the  Laplace  transform  of  this  function  at 
a  discrete  set  of  points  on  the  real  line.  Required  Inputs  are:  N  and  the 
mean  values  of  each  of  these  distributions.  Mean  values  of  these  dist¬ 
ributions  should  be  scaled  so  that  the  sum  of  the  means  Is  near  1. 

Optionally,  a  comparative  Monte-Carlo  simulation  may  be  performed. 

59  'LO'SKIP  2  LINES 

60  PRINT  1  LINE  THUS 
INPUT  THE  VALUE  OF  N. 

62  READ  N 

63  RESERVE  LAV(«)  AS  N 

64  LET  AVGxO.O 

65  LET  VARtO.O 

66  LET  PLAGE: 0 

67  "  IF  FLAGMX  NE  1 

68  "  GO  TO  L3 

69  ' '  OTHERWISE 

70  IF  N  LE  3 

71  LET  FLAGExI 

72  ALWAYS 

73  PRINT  1  LINE  THUS 

INPUT  THE  PROPORTION  OF  THE  1ST  EXPONENTIAL  CONPONENT. 

75  READ  A1 

76  IF  Alxl.O 

77  LET  FLAGMXxO 

78  GO  TO  L3 

79  OTHERWISE 

80  PRINT  1  LINE  THUS 

INPUT  THE  MEAN  VALUE  OF  THE  1ST  EXPONENTIAL  COMPONENT. 

82  READ  THETA 

83  LET  LA lx  1.0/THETA 

84  PRINT  1  LINE  THUS 

INPUT  THE  MEAN  VALUE  OF  THE  2N0  EXPONENTIAL  COMPONENT. 

86  READ  THETA 

87  LET  U2x  1.0/THETA 

88  LET  A2x1.0-A1 

89  LET  LAVCDxLAI 

90  LET  LAV(2)xLA2 

91  LET  AVG1xA1/LAUA2/LA2 

92  LET  VAR1x2.0»(A1/U1»»2*A2/LA2*»2)  -  AVC1«2 

93  LET  AVGiN*AVG1 

94  LET  VARxN4VAR1 

95  GO  TO  L4 

96  'L3'LET  FLACExO 

97  FOR  1x1  TO  N  DO 

98  PRINT  1  LINE  WITH  I 

99  THUS 

INPUT  THE  MEAN  VALUE  OP  THE  ••  TH  RANDOM  VARIABLE. 

D-2 


101  READ  THETA 

102  ADD  THETA  TO  AVG 

103  LET  LAV(I)sCON.AVG/THETA 

104  ADD  CON.AVG»(CON.AVG+1.0)/LAV(I)»»2  -  THETA»»2  TO  VAR 

105  IF  I  >  1 

106  IF  LAV(I)=LAV(I-1) 

107  ADD  1  TO  FLAGE 

108  ALWAYS 

109  ALWAYS 

no  LOOP  ’’OVER  (I)  COMPONENTS 

111  LET  LAULAVd) 

112  IF  FLAGE  s  N-1 

113  LET  FLAGE=1 

114  OTHERWISE 

115  LET  FLAGErO  ’’NOT  ALL  DISTS  ARE  THE  SAME 

116  ALWAYS 

117  *L4*LET  C0ND=1,0/SQRT.F(2.0»PI.C»VAR) 

118  LET  STDV=SQRT,F(VAR) 

119  *•  PRINT  1  LINE  THUS 

120  ”  INPUT  THE  SCALE  FACTOR  (1  GE  GAffU  LE  K2)  IN  LAPLACE  TRANSFORM. 

121  ”  READ  GAFfU 

122  LET  GAFfUsl.O 

123  PRINT  1  LINE  THUS 

DO  YOU  WANT  TO  PERFORM  A  MONTE-CARLO  SIMULATION?  (YES  OR  NO). 

125  READ  ANSWER 

126  IF  SUBSTR.F( ANSWER, 1.1)  *  “Y” 

127  LET  FLAGMrl 

128  PRINT  1  LINE  THUS 

INPUT  THE  INDEX  (1  THRU  9)  OF  THE  RANDOM  #  SEED. 

130  READ  SEED 

131  PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  OF  REPLICATIONS  WANTED. 

133  READ  NREPS 

134  PRINT  1  LINE  WITH  NREPS 

135  THUS 

A  Monte-Cirlo  slaulatlon  of  •••••  roplications  has  bagun. 

137  FOR  Ls1  TO  M,  LET  HISTV(L)*0 

138  LET  AVGT*0.0 

139  LET  VART=0.0 

140  ” 

141  ’’SIMULATE  POR  NREPS  REPLICATIONS. 

142  ” 

143 

144 

145 

146  ” 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 


FOR  1=1  TO  NREPS  DO 
LET  S'jMrO.O 
FOR  Js1  TO  N  DO 

ADD  ERLANG. F(CON.AVG/LAV(J),K, SEED)  TO  SUM 
IF  UNIFORM. F(0.0, 1.0, SEED)  LE  A1 

ADD  EXPONENTIAL. F( 1.0/LAI, SEED)  TO  SUM 
OTHERWISE 

ADD  EXPONENTIAL. F(1.0/LA2, SEED)  TO  SUM 
ALWAYS 

LOOP  ’’OVER  J 
ADD  SUM  TO  AVGT 
ADD  SUM»«2  TO  VART 
FOR  L=  1  TO  M  DO 
IF  SUM  LE  TV(L) 

ADD  1  TO  HISTV(L) 


150 

GO  TO  K2 

159 

OTHERWISE 

160 

LOOP  "OVER  (L)  CELLS 

I6l 

LOOP  "OVER  (I)  REPLICATIONS 

162 

LET  AVOT*AVGT/NREPS 

163 

LET  VART*VART/NREPS-AV0T»»2 

164 

PRINT  1  LINE  THUS 

Monte-Carlo  simulation  has  been  completed. 

166  LET  CDF.MCsO.O 

167  FOR  Lz1  TO  M  DO 

168  LET  PDF«HISTV(L)/NREPS 

169  ADD  PDF  TO  CDF.HC 

170  LET  CDFV(L)sCDF.MC 

171  LOOP  "OVER  (L)  CELLS 

172  OTHERWISE 

173  LET  FLAGMsO 

179  ALWAYS 

175  " 

176  "GET  T'FORHS  PDF  AND  CDF  AT  M  POINTS.  PUCE  IN  DSTARV  t  OSTARV. 

177  •• 

178  FOR  1*1  TO  M  DO 

179  LET  S*CA1#1A*I 

180  LET  CSTARV(I)sLTRNS.FUN  (N,  UV(*),  S) 

181  LET  DSTARV(I)*CSTARV(I)<*S 

182  LOOP  ''OVER  (I)  POINTS  ON  THE  REAL  LINE  IN  THE  S-PUNB 

183  ” 

181  "FILL  THE  ELEMENTS  OF  AM(«,»). 

185  •• 

186  FOR  1*1  TO  M  DO 

187  FOR  J«1  TO  H  DO 

188  LET  E*GAIttt«I  -1.0 

189  LET  AH(I,J)>0.5*WV(J)*(0.5*(XV(J)«1.0))ME 

190  LOOP  "OVER  (J)  COLUMNS 

191  LOOP  "OVER  (I)  ROWS 

192  " 

193  "SOLVE  THE  EQUATIONt  AM  •  GV  .  CSTARV. 

191  " 

195  LET  J>0 

196  CALL  SCBFA  {AM(*,*),  IPVT(«),  J) 

197  IF  J  NE  0 

198  PRINT  1  LINE  WITH  J 

199  THUS 

TROUBLE  FRACTWING  THE  MATRIX  AM  IN  PROGRAM  LP.INV.  J  «  •. 

201  STOP 

202  OTHERWISE 

203  CALL  SGEDI  (AM(*,o),  IPVT(*),  OET(*),  11) 

201  IF  0ET(2)  <  -82 

205  SKIP  2  LINES 

206  PRINT  2  LINES  WITH  0ET(1),  DBT(2) 

207  THUS 

DETERMINANT  OP  MATRIX  AM  «  e.eeee  %  10  EXPOH  (••••),  WHICH  IS  ALMOST 
SINGULAR.  ACCURACY  OF  INV(AH)  IS  OUESTIONABU. 

210  ALWAYS 

211  "  CALL  HAT. INVERSE  (M.  AMI*,*)) 

212  CAU  HAT.VEC.MPT  (AM(*,«).  OSTARVI*},  M)  YIELOINC  CV(«) 

213  CAU  HAT.VEC.MPT  (AM(«.«),  DSTARVI*).  H)  YIELOINC  0V(*) 

211  " 
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y.v.v.v- 


215  "PRINT  OUTPUT  C.D.F. 

216  " 

217  SKIP  2  LINES 

218  PRINT  7  LINES  WITH  N.DIST.NAME 

219  THUS 

(1) 

PROB  DISTRIBUTION  OF  THE  CONVOLUTION  OF  A  SET  OF  «•  DIST'S 


Indep  N-fold  Convolutloii  Norm  Prob  Dlstrlb  Difference 
Variable  p.d.f.  c.d.f.  p.d.f.  c.d.f.  c.d.f. 


227  LET  MAE.DIFF=0.0 

228  LET  HMS.DIFF=0.0 

229  LET  MAE.MCzO.O 

230  LET  RMS.MCrO.O 

231  LET  MAE.NIrO.O 

232  LET  RMS.NI=0.0 

233  FOR  I  BACK  FROM  M  TO  1  DO 

234  SKIP  1  LINE 

235  LET  TrTV(M-I+1) 

236  LET  PDF=DV(I) 

237  LET  CDF=GV(I) 

238  LET  NPDF:C0ND»EXP.F(-0.5*((T-AVG)/STDV)»»2) 

239  LET  NCDF:SMORM((T-AVG)/STDV) 

240  LET  DIFF=CDF-NCDF 

241  LET  MAE.DIFF=riAX.F{MAE.DIFF,ABS.F(DIFF)) 

242  ADD  DIFF»»2  TO  RMS.DIFF 

243  PRINT  1  LINE  WITH  T,PDF,CDF,NPDF,MCDF,DIFF 

244  THUS 

•.••••a*  •.»•••••  t. 

246  IF  FLAGE=1  "ALL  RATE  FARMS  ARE  THE  SAME 

247  "  CALL  NFOLD.U  GIVEN  N,  0.5»LAV(1)»T  TIELDING  EPDF.ECDF 

248  IF  A1  =1.0 

249  CALL  ERLANG  (K«N,  U1,  T)  YIELDING  EPDF.ECDF 

250  OTHERWISE 

251  CALL  NFOLD.MIXE  (N,  A1,  LAI,  LA2,  T)  YIELDING  EPDF.ECDF 

252  ALWAYS 

253  LET  RESID=CDF-ECDF 

254  LET  MAE.NI=MAX.F(MAE.NI,ABS.F(RESID)) 

255  ADD  RESID»»2  TO  BMS.NI 

256  ’Ll'  PRINT  1  LINE  WITH  EPDF.ECDF, CDF -ECDF 

257  THUS 

Exact  fun  ».»»•*»•  •.•»••••  Dlf  rel  to  exact  edf  »,»••*** 


259  OTHERWISE 

260  LET  ECDF=CDF 

261  ALWAYS 

262  LET  DIFF.MC=CDFV(M-I*1)  -  ECDF 

263  LET  MAE.MC=MAX.F(MAE.MC,ABS.F(DIFF.MC)) 

264  ADD  DIFF.MC»»2  TO  RMS.MC 

265  LOOP  "OVER  (I)  CDF  POINTS 

266  PRINT  2  LINES  THUS 


269  LET  RMS.DIFF=SORT.F(BMS.DIFF/REAL.F(M)) 

270  LET  RMS.NI=SQRT.F(RMS.MI/REAL.F(M)) 

271  LET  RMS.MC=SQRT.F(RMS.MC/REAL.F(M)) 
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zn  PRINT  2  LINES  WITH  M,MAE.DIFF,M,HMS.DIFF 

273  THUS 

Max  aba  error  in  Normal  approx  (over  to  c.d.f. 

RMS  difference  In  c.d.f.  and  Normal  approx  (over  ••)  •,•»•••• 

276  PRINT  2  LINES  WITH  MAE.NI,RMS.NI 

277  THUS 

Max  abs  error  in  num  inverse  est  of  the  c.d.f.  •,»•••»• 

RMS  error  in  the  num  inverse  est  of  the  c.d.f. 

280  IF  FLAGM=1 

281  PRINT  2  LINES  WITH  MAE.MC .RMS.MC 

282  THUS 

Max  abs  error  in  Monte-Carlo  estimate  of  the  c.d.f.  •.•••••• 

RMS  error  of  the  Monte-Carlo  estimate  of  the  c.d.f. 

285  ALWAYS 

286  PRINT  3  LINES  WITH  AVO.STDV.DIST.NAME 

287  THUS 

Mean  of  Convolution  Distribution  std  Dev  •»' 

(1)  Mean  of  each  of  the  »••••••••••  distributions: 


291  FOR  Ir1  TO  N  DO 

292  IF  FLAOMXrl 

293  LET  AVGXIrA1/LAV(1)+A2/LAV(2) 

294  OTHERWISE 

295  LET  AVGXI=CON.AVG/LAV(I) 

296  ALWAYS 

297  PRINT  1  LINE  WITH  I.AVGXI 

298  THUS 
Number  ••  Mean 

300  LOOP  "OVER  (I)  COMPONENT  COMPONENTS 

301  SKIP  2  LINES 

302  IF  FLAGM  NE  1 

303  GO  TO  L2 

304  OTHERWISE 

305  PRINT  7  LINES  WITH  N.DIST.NAME.NBEPS 

306  THUS 

SAMPLE  PROS  DIST  OF  THE  SUM  OF  A  SET  OF  ••  •••••••»•»*  RANDOM  VARIABLES 

Monte-Carlo  Sample  »••»» 


Indep 

Variable 

Hlsto  Sample 

Frequency  pedef. 

Sample 
c.d.f • 

LETTcDFiO.O 

315 

FOR  1=1  TO  rt  DO 

316 

LET  XPDF=HISTV(I)/NHEPS 

317 

ADD  XPDF  TO  XCDF 

318 

PRINT  1  LINE  WITH  TV(I) ,HISTV(I) , XPDF, XCDF 

319 

THUS 

331 

LOOP  "OVER  (I)  HISTO  CELLS 

322 

PRINT  2  LINES  THUS 

325  LET  SDT=SQRT.F(VART) 

326  LET  SET=SDT/SQRT.F(REAL.F(NREPS)) 

327  PRINT  3  LINES  WITH  AVGT,SDT,AVGT-1.96»SET,AV0T+1 .96»SET 

328  THUS 


Sample  Average  Value  »»••,»**»  Sample  Standard  Deviation 
95  percent  confidence  interval  in  mean:  •»»»,»»»» 

332  '12* PRINT  1  LINE  THUS 

DO  YOU  HAVE  SIMILAR  PROBLEMS  TO  SOLVE?  (YES  OR  NO). 

334  READ  ANSWER 

335  IF  SUBSTR.F( ANSWER, 1,1)  =  "Y" 

336  RELEASE  LAV(«) 

337  GO  TO  LO 

338  OTHERWISE 

339  STOP 

340  END  "LP.INV 

1  ROUTINE  FOR  MAT.VEC.MPY  (AM,  BV,  NELMTS)  YIELDING  CV 

2  •  * 

3  "ROUTINE  TO  MULTIPLY  THE  SQUARE  MATRIX  AM  ,  OF  NELMTS  BY  NELMTS, 

4  "BY  THE  VECTOR  BV  (NELMTS  BY  1 ) ,  YIELDING  THE  VECTOR  CV  (NELMTS  BY  1 ) . 

5  " 

6  DEFINE  I,  J,  K,  NELMTS  AS  INTEGER  VARIABLES 

7  DEFINE  BV,  CV  AS  REAL,  1-DIMENSIONAL  ARRAYS 

8  DEFINE  AM  AS  A  REAL,  2-DIMENSIONAL  ARRAY 

9  RESERVE  CV(*)  AS  NELMTS 

10  FOR  1=1  TO  NELMTS  DO 

11  LET  CV(I)=0.0 

12  FOR  K=1  TO  NELMTS  DO 

13  ADD  AM(I,K)«BV(K)  TO  CV(I) 

14  LOOP  "OVER  K 

15  LOOP  "OVER  I 

16  RETURN 

17  END  "ROUTINE  MAT.VEC.MPY 

1  FUNCTION  LTRNS.FUN  (N,  LAV,  S) 

2  " 

3  "Obtains  Laplace  transform  of  a  convolution  of  N  prob  dist  functions 

4  "having  rate  parameters  LAV(*).  Complex  argument  (S)  is  evaluated 

5  "only  on  the  real  line.  The  inv  t*form  of  this  function  is  the  c.d.f. 

6  "This  function  must  be  particularized  for  the  desired  form  of  the  prob 

7  "functions.  The  form  used  here  is  indicated  by  the  coounent  statements. 

8  " 

9  DEFINE  I,N  AS  INTEGER  VARIABLES 

10  DEFINE  LAV  AS  A  REAL,  1-DIMENSIONAL  ARRAY 

11  LET  F=1.0/S 

12  " 

13  "CODE  FOR  UNIFORM  WITH  MEAN  =  1/LAV(i). 

14  • » 

15  "  FOR  1  =  1  TO  N,  LET  F=F/S'»(  1 .0-EXP. F(-2.0»S/LAV(I) ) ) 

16  " 

17  "CODE  FOR  EXPONENTIAL. 

18  " 

19  "  FOR  1=1  TO  N,  LET  F=F»LAV(I )/ (S  ♦  LAV(I)) 

20  " 

21  "CODE  FOR  GAMMA(2). 

22  ‘ 

23  "  FOR  1=1  TO  N,  LET  F=F»LAV(I )««2/ (S  ♦  LAV(I))«»2 

24 

25 

26 


CODE  FOR  GAMMA (3). 


-.T  *-1  k-.-i  K^Al.ylL^ 


27  *’  FOR  1=1  TO  N,  LET  F=F«LAV(I)»»3/(S  >  LAV(I))«3 

28  * » 

29  *'CODE  FOR  EXPONENTIAL  MIX. 

30  ** 

31  FOR  1=1  TO  N  DO 

32  IF  A1  =  1.0 

33  LET  F=F»LAV(1)/(S+LAV(1)) 

34  OTHERWISE 

35  LET  F=F« ( A1«LAV(  1 )/ (S^-LAVC  1 ) )+A2»LAV(2)/ (S+LAV(2) ) ) 

36  ALWAYS 

37  LOOP  **OVER  I 

38  RETURN  WITH  F 

39  END  “FUNCTION  LTRNS.FUN 


1  ROUTINE  SGEFA  (  A,  IPVT,  INFO) 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 


“FACTORS  THE  MATRIX  A(»,»)  INTO  UPPER  (U)  AND  STRICTLY  LOWER  (L) 
“TRIANGULAR  MATRICES  SUCH  THAT  A(»,»)  =  U(»,»)L(»,») .  ROUTINE 
' ' IS  INTENDED  FOR  USE  WITH  OTHER  ROUTINES  OF  THE  LINEAR  OPERATIONS 
“PACKAGE— LINPACK.  THIS  VERSION  IS  A  CONVERSION  OF  THE  FORTRAN 

“ROUTINE  WRITTEN  BY  CLEVE  MOLER,  U.  OF  N.M.  AND  ARGONNE  NAT  LAB. 

1 1 

“ARGUMENTS: 

* ' NAME  MODE  ENTRY  VALUE  RETURN  VALUE 

1 1 

'  •  A  HEAL(N,  M)  square  MATRIH  itPPER  TRlAlkdtjLAR  IUTRI]^  AND 

' '  HULTIPLIERS  WHICH  HERE  USED  TO 

"  TO  OBTAIN  IT.  ARE  STORED  IN  L. 

"N  INTEGER  ORDER  OF  THE  MATRIX  A.  DIMENSION  OF  A(*,»). 

"IPVT  INTEGERCN).  VECTOR  OF  PIVOT  INDICES. 

"INFO  INTEGER  INDICATOR.  r  0  FOR  NORMAL  VALUE. 

"  s  K  IF  U(K,K)  EQ  0.0.  THIS 

' '  INDICATES  THAT  SCESL  OR  SGEDI 

"  WILL  DIVIDE  BT  0  IF  CALLED. 


22  DEFINE  I,INF0,J,K,KP1 ,L,N,NM1  AS  INTEGER  VARIABLES 

23  DEFINE  IPVT  AS  AN  INTEGER.  1-DIMENSIONAL  ARRAY 
ZH  DEFINE  A  AS  A  REAL,  2-DIMENSIONAL  ARRAY 

25  ’  ’ 

26  "GAUSSIAN  ELIMINATION  WITH  PARTIAL  PIVOTING. 

27  " 

28  LET  NsDIM.F(IPVT(«)) 

29  LET  INFOrO 

30  LET  NMWN-1 

31  IF  NM1  <  1 

32  GO  TO  L7 

33  OTHERWISE 

3«  FOB  K=1  TO  NMl  DO 

35  LET  KP1=K+1 

36  " 

37  "FIND  L  r  PIVOT  INDEX  IN  THIS  COLUMN. 

38  " 

39  LET  SMAX:ABS.F(A(K.K)) 

90  LET  LsK 

91  FOR  I=K*1  TO  N  DO 

92  IF  ABS.F(A(I,K))  >  SMAX 

93  LET  L=I 
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44  LET  SMAX=ABS.F(A(I,K)) 

45  ALWAYS 

46  LOOP  "FOR  MAX  ELEMENT 

47  LET  IPVT(K)=L 

48  " 

49  "ZERO  PIVOT  IMPLIES  THIS  COLUMN  ALREADY  TRIANGULARIZED. 

50  " 

51  IF  A(L,K)  =  0.0 

52  GO  TO  L4 

53  OTHERWISE 

54  " 

55  "INTERCHANGE  IF  NECESSARY. 

56  " 

57  IF  L  =  K 

58  GO  TO  LI 

59  OTHERWISE 

60  LET  TsA(L,K) 

61  LET  A(L,K):A(K,K) 

62  LET  A(K,K)rT 

63  'Ll*  LET  T=-1.0/A(K,K) 

64  FOR  I=K+1  TO  N,  LET  A(I,K)sT»A(I,K) 

65  " 

66  "ROW  ELIMINATION  WITH  COLUMN  INDEXING. 

67  " 

68  FOR  JrKPI  TO  N  DO 

69  LET  TsA(L,J) 

70  IF  LsK 

71  GO  TO  L2 

72  OTHERWISE 

73  LET  A(L,J):A(K,J) 

74  LET  A(K,J)rt 

75  'LE'  FOR  IsK+l  TO  N,  LET  A(I,J)=T*A(I,K)*A(I,J) 

76  LOOP  "OVER  (J)  COLUMNS 

77  GO  TO  L5 

78  •L4*  LET  INFOrK 

79  'LS'LCOP  "OVER  K 

80  ’LV'LET  IPVT(N)=N 

81  IF  A{N,N)=0.0 

82  LET  INFOsN 

83  ALWAYS 

84  RETURN 

85  END  "SGEFA 

1  ROUTINE  SGEDI  (A,  IPVT,  DET,  JOB) 

2  " 

3  " 

4  "SGEDI  COMPUTES  THE  DETERMINANT  AND  INVERSE  OF  A  MATRIX  USING  THE 

5  "RESULTS  PRODUCED  BY  SGEFA. 

6  " 

7  "ARGUMENTS: 

8  "A(*,*)  THE  REAL  FACTORED  MATRIX  FROM  SGEFA  ON  INPUT.  ON  OUTPUT  THE 

9  "  ~  ARRAY  CONTAINS  THE  MATRIX  INV,  IF  REQUESTED.  ELSE,  UNCHANGED. 

10  "IPVT(«)  THE  INTEGER  PIVOT  VECTOR  FROM  SGEFA. 

11  "JOB  _  AN  INTEGER  SWITCH. 

12  "  ~  X  11  FOR  BOTH  DETERMINANT  AND  INVERSE. 

13  "  X  01  FOR  INVERSE  ONLY. 

14  "  X  10  FOR  DETERMINANT  ONLY. 


15  "DET(»)  CONTAINS  THE  DETERMINANT  OF  THE  MATRIX,  IF  REQUESTED.  ELSE, 

16  "  IS  NOT  REFERENCED.  DETERMINANT  =  DET( 1 )»10.0»»DET{2) ,  WITH 

17  "  DET(I)  BETWEEN  0  AND  10,  AND  WITH  DET(2)  A  FLOATED  INTEGER. 

18  " 

19  "NOTE?  A  DIVISION  BY  ZERO  WILL  OCCUR  IF  THE  INPUT  FACTOR  CONTAINS  A 

20  "ZERO  ON  THE  DIAGONAL  AND  THE  INVERSE  IS  REQUESTED. 

21  " 

22  DEFINE  I,J,J0B,K,KB,KP1,L,N,NM1  AS  INTEGER  VARIABLES 

23  DEFINE  IPVT  AS  AN  INTEGER.  1-DIMENSIONAL  ARRAY 

24  DEFINE  DET,  WORK  AS  REAL,  1-DIMENSIONAL  ARRAYS 

25  DEFINE  A  AS  A  REAL,  2-DIMENSIONAL  ARRAY 

26  LET  N=DIM.F(IPVT{»)) 

27  RESERVE  WORK(»)  AS  N  "LOCALLY 

28  '  • 

29  "CALCULATE  THE  DETERMINANT  IF  REQUESTED. 

30  " 

31  IF  DIV.F(JOB,10)  =  0 

32  GO  TO  L6 

33  OTHERWISE 

34  LET  DETCDsl.O 

35  LET  DET(2)=0.0 

36  LET  TEN* 10.0 

37  FOR  1=1  TO  N  DO 

38  IF  IPVT(I)  NE  I 

39  LET  DET(I)  =  -DET(1) 

40  ALWAYS 

41  LET  DET(1)=A(I,I)«DET(1) 

42  IF  DET(1)=0.0 

43  GO  TO  L6 

44  OTHERWISE 

45  'Ll'  IF  ABS.F(DET(1))  GE  1.0 

46  GO  TO  L2 

47  OTHERWISE 

48  LET  DET(1)sTEN»DET(1) 

49  SUBTRACT  1.0  FROM  DET(2) 

50  GO  TO  LI 

51  'L2'  IF  ABS.F(DET(1))  <  TEN 

52  GO  TO  L4 

53  OTHERWISE 

54  LET  DET(1)*DET(  D/TEN 

55  ADD  1.0  TO  DET(2) 

56  GO  TO  L2 

57  •L4'L00P  "OVER  I 

58  " 

59  "GET  INVERSE  OF  UPPER  TRIANGULAR  MATRIX  U(«,4). 

60  " 

61  •L6'IF  M0D.F(J0B,10)s0 

62  RELEASE  WORK(«) 

63  RETURN 

64  OTHERWISE 

65  FOR  K=1  TO  N  DO 

66  LET  A(K,K):1.0/A(K,K} 

67  LET  Ti-A(K.K) 

68  FOR  1=1  TO  K-  1,  LET  A(I,K)=T*A(I,K) 

69  LET  KP1=K»1 

70  IF  N  <  KP1 

71  GO  TO  L9 


72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 

106 

107 

108 
109 


OTHERWISE 
FOR  JrKPI  TO  N  DO 
LET  T=A(K,J) 

LET  A(K,J)=0.0 

FOR  1=1  TO  K,  LET  A(I,J)=T*A(I,K)+A(I,J) 
LOOP  “OVER  J 

•L9*LOOP  “OVER  K 
( • 

"FORM  INVERSE(U)»INVEHSE{L) 

t  I 

LET  NMUN-1 
IF  NM1  <  1 

RELEASE  WORK(*) 

RETURN 

OTHERWISE 

FOR  KB=1  TO  NM1  DO 
LET  K=N-KB 
LET  KP1=K+1 
FOR  I=KP1  TO  N  DO 

LET  WORK(I)=A(I,K) 

LET  A(I,K)=0.0 
LOOP  “OVER  I 
FOR  JrKPI  TO  N  DO 
LET  TrWORKCJ) 

FOR  1=1  TO  N,  LET  A(I,K)=T»A(I, J)+A(I,K) 
LOOP  • • OVER  J 
LET  L=IPVT(K) 

IF  L  NE  K  “SWAP  ELEMENTS  OF  VECTORS  K  AND  L 
FOR  1=1  TO  N  DO 
LET  T=A(I,K) 

LET  A(I,K)=A(I,L) 

LET  A(I,L)=T 
LOOP  “OVER  I  TO  SWAP 
ALWAYS 

LOOP  “OVER  KB 
RELEASE  WORK(») 

RETURN 
END  “SGEDI 


1  ROUTINE  ERLANG  GIVEN  N,  R,  T  YIELDING  PDF,  CDF 

2  “ 

3  “Calculates  the  probability  density  function  (PDF)  and  cum  dlstrlbu 

4  “tlon  function  (CDF)  for  an  Erlang  function  with  Integer  shape  para 

5  “meter  N,  with  rate  parameter  R,  and  with  real  argument  T. 

6  “ 

7  DEFINE  I,N  AS  INTEGER  VARIABLES 

8  LET  Z=R»T 

9  LET  EXPZ=EXP,F(-Z) 

10  LET  FACTrl.O 

11  LET  ZIrl.O 

12  LETSUMrI.O 

13  FOR  1=1  TO  N-l  DO 

14  LET  FACT=FACT»I 

15  LET  ZI=ZI»Z 

16  ADD  ZI/FACT  TO  SUM 

17  LOOP  “OVER  I 

18  LET  PDF=R»ZI/FACT»EXPZ 
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19  LET  CDF=1.0-EXPZ»SUM 

20  RETURN 

21  END  "ROUTINE  ERLANG 

1  ROUTINE  NFOLD.MIXE  (N,A1  ,LA1  ,U2,T)  YIELDING  PDF, CDF 

2  " 

3  "Routine  produces  the  probability  density  function  (PDF)  and  cum- 

4  "ulatlve  distribution  function  (CDF)  for  an  N-fold  convolution  of  a 

5  "two- component  exponential  mixture  function  having  first  proportion  of 

6  "A1  and  with  rate  parameters  LAI  and  LA2.  Real-valued  argument  Is  T. 

7  " 

8  DEFINE  I,J,K,N  AS  INTEGER  VARIABLES 

9  LET  A2r1.0-A1 

10  LET  EUEXP.F(-LA1»T) 

11  LET  E2=EXP.F(-LA2»T) 

12  IF  Ns1 

13  LET  PDF=a1»LA1»EUA2»LA2»E2 

14  LET  CDF=A1»(1.0-E1)+A2»(1.0-E2) 

15  RETURN 

16  OTHERWISE 

17  IF  Ms2 

18  LET  XCOEFx2.0»A1»A2»LAl»LA2/(LA1-LA2) 

19  LET  PDFx(A1»LA1)»»2»T»EU(A2»LA2)«2«T«E2+XC0EF»(E2-E1) 

20  LET  CDF:1.0-A1»»2»E1»(1.0+LA1«T)-A2«»2»E2»(1.0+LA2«T) 

21  -XCOEF»(E2/LA2-E1/LA1) 

22  RETURN 

23  OTHERWISE 

24  IF  N:3 

25  LET  ARG1:LA1»T 

26  LET  ARG2=LA2»T 

27  LET  F13=1.0-E1»(1.0+AHGU0.5»AR01»*2) 

28  LET  F23=1.0-E2»(1.0*ARG2+0.5«AR02«»2) 

29  "  LET  F12=1.0-E1»(1.0+ARG1) 

30  "  LET  F22=1.0-E2»(1.0+ARG2) 

31  "  LET  FlIrl.O-EI 

32  "  LET  F21=1.0-E2 

33  LET  As1.0/LA1/(LA1-LA2) 

34  LET  B=1.0/LA1/LA2  -  LA1/LA2/(LA2-LA1)»»2 

35  LET  APBcA^B 

36  LET  C=1.0/(LA2-LA1)»»2 

37  LET  AP=1.0/LA2/(LA2-LA1) 

38  LET  BPs1.0/LA1/LA2  -  LA2/LA1/(LA2-LA1)«2 

39  LET  APP:AP+BP 

40  LET  PDF:0.5»A1»»3«LA1»ARG1M2»EU0.5«A2»«3»LA2*ARG2*»2»E2+3.0»A1 

41  •»2»A2»(APB»E1»LA1»»2«LA2-A»LA1«»3*ARG2»EUC»LA1»»2«U2»E2)*3.0 

42  ‘A  1  •A2»»2»  ( APP«E2»U1  •LA2»»2-AP«LA2»«3»ARG  1  •E2+C»LA1  •LA2»»2«E1 ) 

43  LET  CDFrA1«3»F13+A2«3»F23*3.0«A1M2«A2»(1.0-E2»LA1««2/(LA1- 

44  LA2)»»2-(ARG1»LA2/(LA2-LA1)+(U2-2.0«LA1)«LA2/(LA2-LA1)«2)»E1) 

45  ♦ 3 • 0» A 1 • A2»»2» ( 1 . O-E 1 •LA2»»2/ ( LA2-LA 1 )»«2 

46  -(ARC2»LA1/(LA1-LA2)+(U1-2.0»LA2)«LAV(U1-LA2)»»2)«E2) 

47  RETURN 

48  OTHERWISE 

49  PRINT  1  LINE  WITH  N 

50  THUS 

INPUT  ERROR  TO  ROUTINE  NFOLD.HIXE.  N  s 

52  STOP 

53  END  "NFOLD.MIXE 
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ANNEX  E 


SIMSCRIPT  SOURCE  PROGRAM:  INT.TEST 


1  PREAMBLE  '* INT.TEST 

2  NORMALLY  MODE  IS  REAL 

3  DEFINE  SNORM  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

4  DEFINE  ERRFX  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

5  DEFINE  FACTORIAL  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

6  DEFINE  COMPLETE.GAMMA  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

7  DEFINE  NUM.CNVL  AS  A  REAL  FUNCTION  GIVEN  3  ARGUMENTS 

8  END  '* PREAMBLE 

1  MAIN  "INT.TEST 

2  " 

3  "Program  tests  a  variety  of  methods  for  obtaining  convolution  Integral 

4  "of  a  two-parameter  Welbull  distribution.  Program  compares  Leonard 

5  "Johnson’s  approx  for  the  2nd  order  failure  distribution  with  an  exact 

6  "expression  when  time  to  fail  is  a  Welbull  RV.  Ref:  Reliability  and 

7  "Maintainability  of  the  M48A1  Tank,  p.26  ff. 

8  DEFINE  ANSWER  AS  A  TEXT  VARIABLE 

9  DEFINE  FLAGM,I,J,K,KORD,M,N,NCELLS,NREPS,SEED  AS  INTEGER 

10  VARIABLES 

11  DEFINE  HISTV  AS  AN  INTEGER,  1-DIMENSIONAL  ARRAY 

12  DEFINE  TV,FYV,FZV,DELFXV  AS  REAL,  l-DIMENSIONAL  ARRAYS 

13  LET  N=1024  "ELEMENTS  IN  FYV(«) 

14  RESERVE  FYV(«) ,FZV(*) ,DELFXV(»)  AS  N 

15  LET  LINES. V= 9999 

16  LET  NCELLS=20 

17  RESERVE  HISTV (•) ,TV(«)  AS  NCELLS 

18  PRINT  5  LINES  THUS 

Program  calculates  the  convolution  integral  of  N  (N  le  8)  identical  Welbull 
distributions  via  several  methods.  This  convolution  distribution  is  the  c.d.f 
of  the  sum  of  N,  identical  Welbull  random  variables.  Methods  include: 

(a)  evaluation  of  an  analytic  expression,  (b)  Leonard  Johnson's  (L-J)  approxi¬ 
mation,  (c)  finite  numerical  convolution,  and  (d)  Monte-Carlo  simulation. 

24  'LO'SKIP  2  LINES 

25  PRINT  1  LINE  THUS 

INPUT  THE  SCALE  PARAMETER  OF  THE  WEIBULL  DISTRIBUTION. 

27  READ  ETA 

28  LET  Cs2.0/ETA»»2 

29  PRINT  1  LINE  THUS 
INPUT  THE  WEIBULL  SHAPE  PARAMETER. 

31  READ  SHAPE 

32  'Ll 'PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  (LE  8)  OF  CONVOLUTIONS  OF  THIS  DISTRIBUTION  WANTED. 

34  READ  KORD 

35  IF  KORD  >  8 

36  GO  TO  LI 

37  OTHERWISE 

38  LET  ORDER=XORD 

39  ** 

40  "FILL  THE  ARRAYS  OF  DISCRETE  VALUES  OF  THE  C.D.F. 

41  " 

42  LET  FXOrO.O 

43  LET  ERR =0.00001 


U4  LET  TMAX=ETA»(-L0G.E.F(ERH))»»(1. 0/SHAPE) 

45  LET  A VG1=ETA»C0MPLETE,GAtf«( 1.0+ 1.0/SHAPE) 

46  LET  VAR1=ETA»»2»C0MPLETE.GAMMA( 1.0+2. 0/SHAPE)  -  AVG1»»2 

47  LET  STDV1=SQRT.F(VAR1) 

48  LET  TMAX=MAX . F ( TMAX ,ORDER»AVG 1+3 . 5»SQRT . F (ORDERLY AR 1 ) ) 

49  LET  DELZ=TMAX/N 

50  FOR  1=1  TO  N  DO 

51  LET  X=I»DELZ 

52  LET  FXrl.O  -  EXP.F(-(X/ETA)»«SHAPE) 

53  LET  FYV(I)=FX 

54  LET  DELFXV(I)=FX.FX0 

55  LET  FXOsFX 

56  LOOP  »*OVER  (I)  DISCRETE  POINTS  OF  THE  CDF 

57  IF  KORD  >  2 

58  PRINT  1  LINE  WITH  N 

59  THUS 

Starting  numerical  convolution  with  points. 

61  FOR  K=1  TO  K0RD.2  DO 

62  FOR  1=1  TO  N,  LET  FZV(I)=NUM.CNVL  (I,  FYV(»),  DELFXV(»)) 

63  FOR  1=1  TO  N,  LET  FYV(I)=FZV(I) 

64  LOOP  "OVER  (K)  CONVOL  ORDERS 

65  PRINT  1  LINE  THUS 
Numerical  convolution  completed. 

67  ALWAYS 

68  LET  PSI sCOMPLETE .GAMMA (ORDER+ 1 . 0/COMPLETE .GAWA (1.0+1. 0/SHAPE )/ 

69  SHAPE)/FACTORIAL(KORD) 

70  LET  DELT=TMAX/NCELLS 

71  FOR  K=1  TO  NCELLS,  LET  TV(K)=K»DELT 

72  PRINT  1  LINE  THUS 

DO  YOU  WANT  A  MONTE-CARLO  ESTIMATE  OF  THE  CONVOLUTION  C.D.F.?  (Y  OR  N). 

74  READ  ANSWER 

75  IF  SUBSTR.F( ANSWER,  1,1)  =  "Y*» 

76  LET  FLAGM= 1 

77  PRINT  1  LINE  THUS 
INPUT  THE  INDEX  OF  THE  RANDOM  #  SEED. 

79  READ  SEED 

80  PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  OF  REPLICATIONS  WANTED. 


82 

83 

84 


READ  NREPS 

PRINT  1  LINE  WITH  NREPS 
THUS 


A  Monte-Carlo  simulation  of  •••••  replications  has  begun. 


86 

87 

38 

39 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 


FOR  K=1  TO  NCELLS,  LET  HISTV(K)=0 
LET  AVGTsO.O 
LET  VARTrO.O 

I 

»SIfWLATE  FOR  NREPS  REPLICATIONS. 

I 

FOR  1=1  TO  NREPS  DO 
LET  SUMsO.O 
FOR  J=1  TO  KORD  DO 

ADD  W£IBULL.F(SHAP£,ETA»SEED)  TO  SUM 
LOOP  "OVER  J 
ADO  SUM  TO  AVCT 
ADD  SUM«*2  TO  VART 
FOR  K=1  TO  NCELLS  DO 
IF  SUM  LE  TV(K) 


101 

102 

103 

104 

105  *K2* 

106 

107 

108 


ADD  1  TO  HISTV(K) 

GO  TO  K2 
OTHERWISE 

LOOP  **OVER  (K)  CELLS 
LOOP  **OVER  (I)  REPLICATIONS 
LET  AVGT=AVGT/NREPS 
LET  VARTsVART/NREPS.AVGT«»2 
PRINT  1  LINE  THUS 


Monte-Carlo  simulation  has  been  completed. 

110  OTHERWISE 

111  LET  FLAGMrO 

112  ALWAYS 

113  SKIP  2  LINES 

114  PRINT  1  LINE  WITH  KORD, SHAPE, ETA 

115  THUS 

CONVOLUTION  C.D.F.  OF  ORDER  •»  OF  A  WEIBULL  DIST:  SHAPE  AND  SCALE 
117  PRINT  5  LINES  THUS 


Indep  Exact  L-JTpfx  Numerical*  Normal  “His to  Sample 

Variable  c.d.f.(l)  c.d.f.  c.d.f.(2)  Aprx  Freq  c.d.f. 


123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 
157 


LET  CDFOTo 

LET  MAELJ=0.0  ‘TOR  MAX  ABS  ERROR  IN  C.D.F.  FOR  L-J  APPROX 

LET  MAEDNsO.O  '*FOR  MAX  ABS  ERROR  IN  C.D.F  FOB  DISCRETE  NUM  APPROX 

LET  MAEHCrO.O  "FOR  MAX  ABS  ERROR  IN  C.D.F*  FOR  MONTE  CARLO 

LET  MAENA=0.0  "FOR  MAX  ABS  ERROR  IN  C.D.F.  OR  NORMAL  APPROX 

LET  RMSLJsO.O 

LET  RMSDNsO.O 

LET  RMSMCsO.O 

LET  RMSNAsO.O 

LET  AVG.ORDER^AVGI 

LET  VARxORDER*VAR1 

LET  STOVxSQRT.F(VAR) 

FOR  Kx1  TO  NCELLS  DO 
LET  TxTV(K) 

LET  Mx(T40.4999«DELZ)/DELZ 
LET  Xx(PSI»T/ETA) ••SHAPE 
LET  SUMxl.O 
LET  FACTS  1.0 
LET  XXrl.O 

FOR  1=1  TO  KORD-1  DO 
LET  FACTsFACT*! 

LET  XI=XI«X 
ADD  XI /FACT  TO  SUM 
LOOP  "OVER  I 

let  QK=1.0  -  EXP.F{-X)^SUM 
IF  KORDsI 

LET  FZ=FYV(M) 

OTHERWISE 

LET  FZsNUM.CNVL  (M,  FYV(«),  DELFXV(«)) 

ALWAYS 

IF  SHAPEs2.0  AND  K0RD=2 
LET  ARGrT»SQRT.F(C/2.0) 

LET  INTC=EXP.F(0.5»ARC»^2)«(SN0RM(ARG)-SN0RM(-ARG))^ARG^ 
SQRT.F(PI.C/2.0) 

LET  Q2=1.0  -  EXP.F(-ARG«»2)®{1.0*INTC) 
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158  OTHERWISE 

159  LET  Q2=FZ 

160  ALWAYS 

161  LET  MAELJ:MAX.F(MAELJ,ABS.F(Q2-QK)) 

162  LET  MAEDN:MAX.F(MAEDN,ABS.F(Q2-FZ)) 

163  LET  FN=SNORM((T-AVG)/STDV) 

164  LET  NERRsFM-Q2 

165  LET  MAENA:;MAX.F(MAENA,ABS.F(NERR)) 

166  ADD  NERR**2  TO  RMSNA 

167  ADD  (Q2-QK)»»2  TO  RMSLJ 

168  ADD  (Q2-FZ)»»2  TO  RMSDN 

169  IF  FLAOM=1 

170  LET  PDFX=HISTV(K)/NREPS 

171  ADD  PDFX  TO  CDFX 

172  LET  MAEMC=MAX.F(MAEMC,ABS.F(Q2-CDFX)) 

173  ADD  (Q2-CDFX)»»2  TO  RHSMC 

174  PRINT  1  LINE  WITH  T,Q2,QK,FZ,FN,HISTV(K),CDFX 

175  THUS 

!.•••••■ 

177  OTHERWISE 

178  PRINT  1  LINE  WITH  T,Q2,QK,FZ,FN 

179  THUS 

181  ALWAYS 

182  LOOP  ‘'OVER  (K)  VALUES  OF  TIME 

183  LET  RMSLJ:SQRT.F(RHSLJ/REAL.F(NCELLS)) 

184  LET  RMSON:SQRT.F(RMSOH/REAL.F(NCELLS)) 

185  LET  RMSHC=SQRT.F(RMS^«:/REAL.F(NCELLS)) 

186  LET  RMSNA<SQRT.F(RMSNA/REAL.F(NCELLS)) 

187  PRINT  4  LINES  WITH  M 

188  THUS 

rn  fhe~drscr«t*  "nmifiSi~Spff6x  exact  it  tiilMr  the 

WelbuU  shape  parameter  Is  not  2  or  the  number  of  convolutions  Is  not  2. 
(2)  Number  of  discrete  points  In  numerical  convolution 

193  PRINT  4  LINES  WITH  MAELJ.RMSU.MAEON.RNSDN.MAENA.RHSNA.HAEMC, RHSMC 

194  THUS 

Max  abs  error  and  RMS  error  In  c.d.f.  of  L*J  approximation  e.eeeeee  a,asaaae 

Max  abs  error  and  RMS  error  In  c.d.f.  of  discrete  numerical  e.eeeeee  e.eeeeee 

Max  abs  error  and  RMS  error  In  c.d.f.  of  Normal  approx  a.aeeeae  e.eeaaaa 

Max  abs  error  and  RMS  error  In  c.d.f.  of  Monte>Carlo  aim  e.eeeeee  e.eeeeee 

199  PRINT  2  LINES  WITH  AVO.STDV 

200  THUS 

Mean  of  ine  convolution  distribution  eeee.eeee  std  Dev  eeee.eeee 

203  IF  FLAGMrl 

204  let  SET=SQRT.F(VART/REAL.F(MREPS)) 

205  PRINT  3  LINES  WITH  AWCT,SORT.r(VART),AVOT-1.96»SET,AVCT*1.96»SBT 

206  THUS 

Sample  average  of  sum  of  Welbull  RVs  sees. sees  gtd  Dev  eeee.eeee 

95  percent  confidence  Interval  In  mean  eeee.eeee.  eeee.eeee 

210  ALWAYS 

211  PRINT  1  LINE  THUS 

DO  YOU  WANT  TO  CONTINUE?  (YES  OR  HO). 

213  READ  ANSWER 
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214 

IF  SUBSTR,F( ANSWER, 1,1)  s 

tiyn 

215 

GO  TO  LO 

216 

OTHERWISE 

217 

STOP 

218 

END  "INT.TEST 

1 

o 

FUNCTION  FACTORIAL(N) 

1 1 

c 

3 

"CALCULATES  THE  FACTORIAL  OF 

INTEGER 

4 

« 1 

5 

DEFINE  I  AND  N  AS  INTEGER  VARIABLES 

6 

IF  N  LE  0 

7 

RETURN  WITH  1.0 

8 

OTHERWISE 

9 

LET  FrI.O 

10 

FOR  Ir1  TO  N,  LET  F=F»I 

11 

RETURN  WITH  F 

12 

END  "FACTORIAL 

1  FUNCTION  COMPLETE. GA»U (XX) 

2  " 

3  "CALCULATES  THE  COMPLETE  Gxmk  FUNCTION  VIITH  SINGLE  REAL  ARGUMENT  XX. 
<4  "  METHOD!  THE  RECURSION  RELATION  AND  POLXNOMIAL  APPROXIMATION  IS  TAKEN 

5  "FROM:  C.  HASTINGS,  JR,  'APPROXIMATIONS  FOR  DIGITAL  COMPUTERS,' 

6  "PRINCETON  UNIV.  PRESS,  1953. 

7  " 

8  IF  XX  >  57.0 

9  GO  TO  L130 


10 

OTHERWISE 

11 

*L6* 

LET  X  s  XX 

12 

LET  ERR  :  0.000001 

13 

LET  GA^  s  1.0 

14 

IF  X  LE  2.0 

15 

GO  TO  L50 

16 

OTHERWISE 

17 

GO  TO  L15 

18 

*L10* 

IF  X  LE  2.0 

19 

GO  TO  L110 

20 

OTHERWISE 

21 

*L15* 

SUBTRACT  1.0  FROM  X 

22 

LET  Gkmk  3  GAMHA  •  X 

23 

GO  TO  L10 

24 

'L50' 

IF  X  s  1.0 

25 

CO  TO  L120 

26 

OTHERWISE 

27 

IF  X  >  1.0 

28 

GO  TO  L110 

29 

OTHERWISE 

30 

*L60' 

IF  X  >  ERR 

31 

GO  TO  L80 

32 

OTHERWISE 

33 

UT  Y  «  REAL.F(TRUNC.F(X))*X 

34 

IF  ABS.F(Y)  LE  ERR 

35 

GO  TO  L130 

36 

OTHERWISE 

37 

IF  Y^ERR  GE  1.0 

38 

GO  TO  L130 

39 

OTHERWISE 

40 

'L70' 

IF  X  >  1.0 

41 

GO  TO  L110 

42 

OTHERWISE 

43 

•L80' 

LET  GAWU  :  aHMk  /  X 

44 

ADD  1.0  TO  X 

45 

GO  TO  L70 

46 

'L110' 

LET  Y  s  X  -  1.0 

47 

LET  GY  =  1.0*Y»(-0.57710l66*Y»(0.98585399*Y»(-0.876lt2ie2+Y» 

48 

(0.83282120*Y»(-0. 56847290+ Y»(0.25482049»Y»(-0. 05  H*9930))))))) 

49 

LET  GAM1A  ::  GAWU  •  GY 

50 

•L120' 

RETURN  WITH  GAWU 

51 

'L130' 

PRINT  1  LINE  WITH  XX  THUS 

ERROR 

IN  COMPLETE. GAffiA.  ARGUMENT  = 

53 

STOP 

54 

END  " 

COMPLETE.  GAWU 

1  FUNCTION  NUN.CNVL  (N,  FYV,  DELFXV) 

2  ** 

3  ''Function  calculates  a  value  of  the  c.d.f.  of  the  sua  of  two  randoa 

4  "variables—  x  and  y^-whose  c.d.f.'s  are  evaluated  at  a  dlaorete  # 

5  "of  points  on  their  domains.  This  distribution  of  the  sum  ia  the  con* 

6  "volution  of  the  dist's  of  x  and  y.  The  convolution  distribution  is 

7  "evaluated  for  the  N  th  discrete  argument.  The  set  of  c.d.f.  values 

8  "of  y  are  given  by  the  vector  FYV»  and  the  first  backward  differences 

9  "in  the  c.d.f.  of  x,  defined  on  the  same  finite  domain,  are  given  in 

10  "DELFXV. 

11  DEFINE  I.N  AS  INTEGER  VARIABLES 

12  DEFINE  FYV, DELFXV  AS  REAL,  UDIHENSIONAL  ARRAYS 

13  LET  GN»0.0 

14  FOR  Is1  TO  N-l,  ADD  FYV(N-I )»DBLFXV(I)  TO  GN 

15  RETURN  WITH  GN 

16  END  "FUNCTION  NUM.CNVL 


ANNEX  F 


SIMSCRIPT  SOURCE  PROGRAM:  TEST.CONVOLV 


2  PREAMBLE  ’’TEST.CONVOLV 

3  NORMALLY  MODE  IS  REAL 

4  DEFINE  SNORM  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

5  DEFINE  ERRFX  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

6  DEFINE  COMPLETE. GAMMA  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

7  DEFINE  W2WFUN  AS  A  REAL  FUNCTION  GIVEN  1  ARGUMENT 

8  DEFINE  EFUN  AS  A  REAL  FUNCTION  GIVEN  2  ARGUMENTS 

9  DEFINE  WFUN  AS  A  REAL  FUNCTION  GIVEN  2  ARGUMENTS 

10  END  ’’PREAMBLE 

1  MAIN  ’’TEST.CONVOLV 

2  ” 

3  ’’Program  to  run  routine  CONVOLV.  This  program  generates  2  p.d.f.’s 

4  ’’defined  on  a  discrete  point  set,  from  prob  dist’s  to  be  numerically 

5  ’’convolved  via  Fourier  transformation,  multiplication  of  transforms, 

6  ’’and  inversion.  The  #  of  real  (as  opposed  to  imaginary)  points  in 

7  ’’the  transform  and  its  inverse  must  be  a  power  of  2  in  order  to  use 

8  ’’the  Cooley-Tukey  FFT  algorithm.  Comparisons  with  exact  results  and, 

9  ’’optionally,  Monte-Carlo  results  are  also  given, 

10  ” 

DEFINE  FLAGE,FLAGM,FLAGW,I,J,K,L,M,MINCR,N,NFOLD,NCELLS,NREPS,SEED  AS 

12  INTEGER  VARIABLES 

13  DEFINE  ANSWER, FUN. NAME  AS  TEXT  VARIABLES 

14  DEFINE  HISTV  AS  AN  INTEGER,  1-DIMENSIONAL  ARRAY 

15  DEFINE  TV,XV,YV,PDFV,CDFV  AS  REAL,  l-DIMENSIONAL  ARRAYS 

16  PRINT  11  LINES  THUS 

This  program  calculates  the  probability  distribution  of  the  sum  of  a 
set  of  random  variables  of  a  particular  type,  such  as  Erlang  or  Welbull. 

This  is  equivalent  to  obtaining  the  N-fold  convolution  of  the  probability 
functions  of  the  set  of  N.  For  a  given  type  of  random  variable,  two  sets 
of  parameters  are  permitted.  Distributions  having  the  1st  parameter  set  are 
convolved  N-1  times  with  the  distribution  having  the  2nd  parameter  set. 

Where  available,  exact  results  are  calculated  and  displayed.  A  numerical 
method  for  obtaining  convolution  integrals  based  on  the  Fourier  transform 
is  used  in  all  cases  to  obtain  an  approximation  of  Lhe  convolution  p.d.f. 
and  c.d.f.  Optionally,  Monte-Carlo  simulation  is  used  for  sample  estimates. 


28 

PRINT  3  LINES 

29 

THUS 

The  current  program  version  treats  convolutions  of  an  Erlang  or  a 

Weibull 

distribution  in  standardized  form,  i.e.,  characterized  by  a  shape 

parameter. 

33 

PRINT  1  LINE  THUS 

IF  THE 

ERLANG  FORM  IS  WANTED,  INPUT  AN  E;  OTHERWISE,  INPUT  A 

w. 

35 

READ  ANSWER 

36 

IF  SUBSTR.F(ANSWER,1,1)  =  ”E” 

37 

LET  FLAGW=0 

38 

LET  FLAGE=1  ’’TRIGGER  FORMAT  FOR  EXACT  RESULTS 

39 

LET  FUN.NAMEz  ’’Erlang” 

40 

OTHERWISE 

41 

LET  FLAGWrI 

‘*2  LET  FUM.NAMEr  "Weibull" 

43  ALWAYS 

44  PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  OF  CONVOLUTIONS  WANTED. 

46  HEAD  NFOLD 

47  "  LET  N=4096 

48  LET  N=8192 

49  LET  M=DIV.F(N,2)  "NUMBER  OF  REAL  POINTS  IN  THE  SERIES 

50  RESERVE  XV(»),  YV(«)  AS  N 

51  LET  MINCR  =  256  "SKIP  INTERVAL  FOR  PRINTING 

52  LET  NCELLS=DIV.F(M,MINCR) 

53  LET  NCELLS=MAX.F(16,NCELLS) 

54  LET  MINCR=DIV.F(M,NCELLS) 

55  RESERVE  HISTVC*)  AS  NCELLS 

56  RESERVE  TV(»)  AS  NCELLS  "FOR  INDEPENT  VAR  IN  A  HISTOGRAM 

57  RESERVE  PDFVC*) ,CDFV(»)  AS  NCELLS 

53  LET  LINES. V:9999 

59  LET  ETAlrl.O 

50  LET  ETA2=1.0 

61  IF  FLAGW=1 

62  GO  TO  L7 

63  OTHERWISE 

64  'LO’PRINI  1  LINE  WITH  FUN. NAME 

65  THUS 

INPUT  THE  INTEGER  SCALE  PARAM  OF  THE  1ST  STD  »»»»»»»  DISTRIBUTION. 

67  READ  K 

68  IF  K  <  1 

69  PRINT  1  LINE  THUS 

Try  again  using  a  positive  integer. 

71  GO  TO  LO 

72  OTHERWISE 

73  " 

74  "CALCULATE  MEAN  AND  VARIANCE  OF  1ST  DIST. 

75  " 

76  LET  AVGUK 

77  LET  VAR1=K 

78  'Ll’ PRINT  1  LINE  WITH  FUN. NAME 

79  THUS 

INPUT  THE  INTEGER  SCALE  PARAM  OF  THE  2ND  STD  •»»»»»»  DISTRIBUTION. 

81  READ  L 

82  IF  L  <  1 

83  PRINT  1  LINE  THUS 

Try  again  using  a  positive  integer. 

85  00  TO  LI 

36  OTHERWISE 

37  ' ' 

38  "CALCULATE  MEAN  AND  VAR  OF  2ND  DIST  AND  OF  CONVOLUTION  DIST. 

89  ' ' 

90  LET  AVG2:L 

91  LET  VAR2=L 

92  LET  AVU=AV01»(NF0LD-1)+AVG2 

93  LET  VAR=VAH1»(NF0LD-1)+VAR2 

94  LET  STDV=SORT.F(VAR) 

95  LET  STDV1=SURT.F(VAR1) 

96  LET  STDV2=S'JRT.F(VAR2) 

97  SKIP  2  LINES 

98  PRINT  7  LINES  WITH  FUN. NAME, K,L,M 


99  THUS 

EXACT  CONVOLUTION  OF  TWO  *******  DENSITIES  WITH  SHAPE  PARAMS  **  AND  ** 

Number  of  real  points  in  the  Fourier  transform  **** 

Indep  1st  Ist  STd  2nd  ^nv  Conv  Nofmaf 

Variable  p.d.f.  c.d.f.  p.d.f.  c.d.f.  p.d.f.  c.d.f.  c.d.f. 

107  LET  RANGE=AVG+3r7»SQRT.TrVARr 

108  'LU'LET  DARG=RANGE/M 

109  ' ' 

110  "CHECK  RANGE  AND  MODIFY,  AS  NECESSARY. 

111  " 

112  IF  EFUN((NF0LD-1)»K+L, RANGE)  <  0.9999 

113  ADD  DARG  TO  RANGE 

114  GO  TO  L4 

115  OTHERWISE 

116  GO  TO  L8 

117  " 

118  "GET  INPUTS  FOR  WEIBULL  DISTRIBUTION. 

119  " 

120  'L7' PRINT  1  LINE  WITH  FUN. NAME 

121  THUS 

INPUT  THE  SHAPE  PARAMETER  OF  THE  1ST  »»»»»»»  DISTRIBUTION. 

123  READ  SHAPE1 

124  PRINT  1  LINE  WITH  FUN. NAME 

125  THUS 

INPUT  THE  SHAPE  PARAMETER  OF  THE  2ND  »»»»»»»  DISTRIBUTION. 

127  READ  SHAPE2 

128  LET  ERRrO. 00001 

129  LET  T1MAX:ETA1»(-L0G.E.F(ERR))»»{1.0/SHAPE1) 

130  LET  T2MAXiETA2»(-LOG.E.F(ERR))»»(1.0/SHAPE2) 

131  LET  AVGUETA1»C0MPLETE.GAMMA(1.0+1.0/SHAPE1) 

132  LET  AVG2=ETA2»COMPLETE.GAMMA(1.0+1.0/SHAPE2) 

133  LET  VAR1rETA1»»2»COMPLETE.GAMMA(1.0+2.0/SHAPE1)  -  AVG1»»2 

134  LET  STDV1=SQRT.F(VAR1) 

135  LET  VAR2=ETA2»»2»C0MPLETE. GAMMA (1. 0+2. 0/SHAPE2)  -  AVG2»»2 

136  LET  STDV2=SQRT.F(VAR2) 

137  LET  AVG=(NF0LD-1)»AVG1+AVG2 

138  LET  VARr(NFOLD-1)»VAR1+VAR2 

139  LET  STDV=SQRT.F(VAR) 

140  LET  RANGE=MAX.F(T1MAX,T2MAX) 

141  LET  RAN0E=MAX.F(RAN0E,AV0+3.7»SQRT.F(VAR)) 

142  ' ' 

143  "PRINT  HEADINGS  FOR  INPUT  DISTRIBUTIONS. 

144  ' ' 

145  SKIP  2  LINES 

146  PRINT  3  LINES  WITH  FUN. NAME, SHAPE1 ,SHAPE2,M 

147  THUS 

EXACT  CONVOLUTION  OF  TWO  »»•»»»»  DENSITIES  WITH  SHAPE  PARAMS  »».»  AND  »».» 
Number  of  real  points  in  the  Fourier  transform  »••• 


151 

IF  SHAPE1=2.0  AND 

SHAPE2=2.0 

152 

IF  ETAl  NE  ETA2  OR  NFOLD 

153 

GO  TO  L9 

154 

OTHERWISE 

155 

LET  FLAGErl 

156 

PRINT  4  LINES 

THUS 

F-3 


Indep”  fsT"  fsE  2n3  ?n3’  Conv  '  Conv  Normal 

Variable  p.d.f.  c.d.f.  p.d.f.  o.d.f.  p.d.f.  c.d.f.  c.d.f. 

“Tei  goYoTs 

162  OTHERWISE 

163  ’L9’LET  FLAGE=0 

164  PRINT  4  LINES  THUS 


Indep  fst  Ist  2nd  ”  2nd 

Variable  p.d.f.  c.d.f.  p.d.f.  c.d.f. 


■"69  ■ 

■•Iff’ LET 

ARGO=0.0  ’’FIXED 

170 

LET 

RANGE=RANGE-ARGO 

171 

LET 

DARG=RANGE/M 

172 

LET 

AVGTDrO.O  ’’THEORETICAL,  DISCRETIZED 

173 

LET 

VARTDrO.O 

174 

LET 

XSUMr 1 . 0 

175 

LET 

XXSUM=0.0 

176 

LET 

FlOrO.O 

177 

LET 

F20=0.0 

178 

LET 

F30=0.0 

179 

LET 

J  =  0  ”T0  COUNT  CELLS 

180 

LET 

MAENArO.O 

181 

LET 

RMSNArO.O 

182 

1 1 

183 

••GET  TEST  FUNCTIONS. 

184 

1 1 

185 

FOR 

1=1  TO  M  DO 

186 

LET  ARG=I*»DARG>ARGO 

187 

IF  MOD.F(I,2)=0 

188 

LET  COEF=2.0 

189 

OTHERWISE 

190 

LET  C0EF=4.0 

191 

ALWAYS 

192 

IF  FLAGWrO 

193 

LET  F1=EFUN(K,AR0) 

194 

LET  F2=EFUN(L,ARG) 

195 

LET  F3:EFUN((NF0LD-1)»K+L,ARG) 

196 

OTHERWISE 

197 

LET  F1:WFUN(SHAPE1,ARG) 

198 

LET  F2=WFUN(SHAPE2,ARG) 

199 

IF  FLAGE=1 

200 

LET  F3=W2WFUN(ARG) 

201 

ALWAYS 

202 

ALWAYS 

203 

LET  XV(2»I-1)=F1-F10 

204 

LET  YV(2»I-1)=F2-F20 

205 

LET  CDENS=F3-F30 

206 

ADD  CDENS^ARG  TO  AVGTD 

207 

ADD  CDENS*»ARG«»2  TO  VARTD 

208 

LET  FIOrFI 

209 

LET  F20=F2 

210 

LET  F30=F3 

211 

LET  UPPER=1.0-F3 

212 

ADD  COEF<*UPPER  TO  XSUM 

213 

ADD  COEF»ARG»UPPER  TO  XXSUM 

F-4 

214  • 

• 

215  ’ 

•FILL  IMAGINARY  COMPONENTS  WITH  ZEROS. 

216  • 

1 

217 

LET 

XV(2»I)=0.0 

218 

LET 

YV(2»I)=0.0 

219 

IF  1 

M0D.F(I,MINCR)=0 

220 

ADD  1  TO  J 

221 

LET  TV{J)=ARG 

222 

IF  FLAGE=1 

223 

LET  PDFV(J)=CDENS 

224 

LET  CDFV(J)=F3 

225 

LET  FN=SNORM((ARG-AVG)/STDV) 

226 

LET  NERRrFN-F3 

227 

LET  MAENA=MAX . F (MAENA , ABS . F ( NERR ) ) 

228 

ADD  NERR»»2  TO  RMSNA 

229 

PRINT  1  LINE  WITH  ARC ,XV(2«I-1 ) ,F1 

230 

F3,FN  THUS 

........  ........  ........  ........ 

232 

OTHERWISE 

233 

PRINT  1  LINE  WITH  ARG,XV(2»I- 1 ) ,F1 

234 

THUS 

tt.tttVttttll  tt 


««««»«« 


236  ALWAYS 

237  ALWAYS 

238  LOOP  "OVER  (I)  DATA  POINTS 

239  PRINT  2  LINES  THUS 


242  PRINT  4  LINES  WITH  FUN. NAME. AVG1 ,STDV1 , FUN. NAME, AVG2,STDV2.AVG,STDV 

243  THUS 

Mean  and  standard  deviation  of  the  1st  »»»»»»»  dist*ni  »»•,»»»» 

Mean  and  standard  deviation  of  the  2nd  »»»»»»»  dist’n:  •««.«««« 

Theoretical  mean  and  SD  of  the  convolution  distribution;  •»»,»»»»  »••,»»»• 

248  IF  FLAGErl 

249  LET  VARTDrVARTD-AVGTD»»2 

250  LET  XSUM=XSUM»DARG/3.0 

251  LET  XXSUM=2.0»DARG/3.0»XXSUM  -  XSUM»»2 

252  PRINT  3  LINES  WITH  AVGTD, SORT. F ( VARTD ) ,XSUM, SORT. F(XXSUM) 

253  THUS 

Avg  and  SO  of  theoretical,  discretized  convolution  dist;  »»».»»»•  •■».»»»• 

Alternate  (2nd  order)  calculation  of  average  and  std  dev;  »»»,»»»•  »■»,§»»» 

257  LET  RMSNA=SQRT.F(RMSNA/REAL.F(NCELLS)) 

258  PRINT  3  LINES  WITH  MAENA.RMSNA 

259  THUS 

Max  abs  error  and  RMS  error  in  c.d.f.  of  Normal  approx;  ».»»»»»» 


263  ALWAYS 

264  •• 

265  “TAKE  NUMERICAL  CONVOLUTION. 

266  •  • 

267  PRINT  2  LINES  WITH  NFOLD,M 

268  THUS 

Starting  ••  convolutions  using  FFT  with  »»»»  real  points. 


F-5 


V  .  ■ 


271 

CALL  CONVOLV  (NFOLD,  XV(«), 

YV(*»)) 

272 

SKIP  1  LINE 

273 

LET  SUMz  0.0 

274 

FOR  1=1  TO  M,  ADD  YV(2»I-1) 

TO  SUM 

275 

PRINT  1  LINE  WITH  2.0»SUM/N 

276 

THUS 

Cumulative  of  numerical  convolution  density  Is  . 

278  SKIP  2  LINES 

279  PRINT  6  LINES  WITH  NFOLD, FUN. NAME 

280  THUS 

EXACT  VERSUS  NUMERICAL  CONVOLUTION  OF  »»  •••••»»  PROB  DISTRIBUTIONS 


Indep  Theory  Theory  Numer  Numer  Dlff 

Variable  p.d.f.  c.d.f.  p.d.f,  c.d.f.  c.d.f. 

~sf  ^LET'^rffr'oro 

288  LET  J=0 

289  LET  AVGFT=0.0 

290  LET  VARFTrO.O 

291  LET  MAEFTzO.O 

292  LET  RMSFTrO.O 

293  LET  XSUMrI.O 

294  LET  XXSUMrO.O 

295  FOR  1=1  TO  M  DO 

296  LET  ARG=I»DARG4.ARG0 

297  IF  M0D.F(I,2)=0 

293  LET  C0EF=2.0 

299  OTHERWISE 

300  LET  C0EF=4.0 

301  ALWAYS 

302  LET  P0F.FT=YV(2»I-1)»2.0/N 
30 ADD  ARG-PDF.FT  TO  AVGFT 

304  ADD  ARG*»2»PDF.FT  TO  VARFT 

305  ADD  PDF. FT  TO  CDF. FT 

30b  LET  UPPERr1.0-CDF.FT 

307  ADD  COEF»UPPER  TO  XSUM 

308  ADD  COEF»ARG»UPPER  TO  XXSUM 

309  IF  MOD.F(I,MINCR)=0 

310  ADD  1  TO  J 

311  IF  FLAGErl 

312  LET  PDF=PDFV(J) 

313  LET  CDF=CDFV(J) 

314  OTHERWISE 

315  LET  PDFrPDF.FT 

316  LET  CDFrCDF.FT 

317  LET  PDFV(J)=PDF 

318  LET  CDFV(J)rCDF 

319  ALWAYS 

320  LET  DIFFrCDF.FT-CDF 

321  LET  MAEFT=MAX.F(MAEFT,ABS.F(DIFF)) 

322  ADD  DIFF»»2  TO  RMSFT 

323  PRINT  1  LINE  WITH  ABG,PDF,CDF,PDF.FT,CDF.FT,DIFF 

324  thus 

■i.ttii  ■.■•■■■■  I. 

326  ALWAYS 

327  LOOP  ‘’OVER  (I)  POINTS 


328 


PRINT  2  LINES  THUS 


331  LET  VARFT=VARFT-AVGFT»«2 

332  LET  XSUM:XSUM»DARG/3.0 

333  LET  XXSUM=2.0»DARG/3.0»XXSUM  -  XSUM«»2 

334  LET  RMSFT=SQRT.F(RMSFT/REAL.F(NCELLS)) 

335  PRINT  5  LINES  WITH  NFOLD, FUN. NAME, AVGFT, SORT. F(VARFT) ,XSUK, 

336  SQRT.F(XXSUM),MAEFT,RMSFT 

337  THUS 

Mean  value  and  standard  deviation  of  the  suu  of  ••  •»»»»»»  rvs  via  FFT: 

Calculated  mean  »»».••»»  Std  Dev  **»,»«»» 

Alternate  mean  ••».»»»»  std  Dev  •»».»••• 

Max  abs  error  and  RMS  error  In  oonvol  c.d.f.  via  FFT  ».»»»»»»  ».»»»»•» 


343  PRINT  1  LINE  THUS 

DO  YOU  WANT  TO  PERFORM  A  MONTE-CARLO  SIMULATION?  (YES  OR  NO). 

345  READ  ANSWER 

346  IF  SUBSTR.F(ANSWER,1,1)  NE  "Y” 

347  GO  TO  L5 

348  OTHERWISE 

349  "  LET  FLAGMzl 

350  PRINT  1  LINE  THUS 

INPUT  THE  INDEX  (1  THRU  9)  OF  THE  RANDOM  #  SEED. 

352  READ  SEED 

353  PRINT  1  LINE  THUS 

INPUT  THE  NUMBER  OF  REPLICATIONS  WANTED. 

355  READ  NREPS 

356  PRINT  1  LINE  WITH  NREPS 

357  THUS 

A  Monte-Carlo  simulation  of  replications  has  begun. 

359  LET  AVGT=0.0 

360  LET  VARTsO.O 

361  FOR  1=1  TO  NCELLS,  LET  HISTV<I)=0 

362  ” 

363  "SIMULATE  FOR  NREPS  REPLICATIONS. 


FOR  1=1  TO  NREPS  DO 
LET  S0M=0.0 
FOR  J=1  TO  NFOLD-1  DO 
IF  FLAGW=0 

ADD  ERLANG. F(AVG1,K, SEED)  TO  SUM 
OTHERWISE 

ADD  WEIBULL.F(SHAPE1,ETA1,SEED)  TO  SUM 
ALWAYS 

LOOP  "OVER  (J)  RV’S 
IF  FLAGWrO 

ADD  ERLANG. F(AVG2,L, SEED)  TO  SUM 
OTHERWISE 

ADD  WEIBULL.F(SHAPE2,ETA2,SEE0)  TO  SUM 
ALWAYS 

ADD  SUM  TO  AVGT 
ADD  SUM»»2  TO  VART 
FOR  J=1  TO  NCELLS  DO 
IF  SOM  LE  TV(J) 

ADD  1  TO  HISTV(J) 

GO  TO  K2 


385  OTHERWISE 

386  LOOP  "OVER  (J)  CELLS 

387  'KP'LOOP  "OVER  (I)  REPLICATIONS 

388  LET  AVGT=AVGT/NREPS 

389  LET  VART=VART/NREPS  -  AVGT«»2 

390  LET  SET=SQRT.F(VART/REAL.F(NREPS)) 

391  PRINT  3  LINES  WITH  NREPS, AVGT,SQRT.F(VART) ,AVGT-1 .96»SET, AVGT+1 .96»SET 

392  THUS 

Sample  mean  from  »»»»  reps  of  Monte-Carlo  slm  »«».»»»»  std  Dev  »»».*»»• 

93t  statistical  confidence  interval  in  mean:  •»».»»»»,  »»«,»»»» 

396  SKIP  2  LINES 

397  PRINT  6  LINES  WITH  NFOLD,  FUN. NAME 

398  THUS 

SAMPLE  PROB  DIST  OF  THE  SUM  OF  A  SET  OF  »»  »»»»•»»  RANDOM  VARIABLES 


fndep  Hfsto  Sample  Sampfe  Theory  Differ 

Variable  Frequency  p.d.f.  c.d.f.  c.d.f.  c.d.f. 


«05  LET  XCDF=0.0 

U06  LET  MAEMCrO.O 

907  LET  RMSMC=0.0 

903  FOR  J:1  TO  NCELLS  DO 

909  LET  XPDF=HISTV(J)/NREPS 

910  ADD  XPDF  TO  XCDF 

911  LF.T  CDF=CDFV(J) 

912  LET  DIFF=XCDF-CDF 

413  LET  MAEMC=MAX,F(MAEMC,ABS.F(DIFf)) 

414  ADD  DIFF»»2  TO  RMSMC 

415  PRINT  1  LINE  WITH  TV(J) ,HISTV(J) , XPDF, XCDF, CDF, DIFF 

416  THUS 

«««•»» 

418  LOOP  "OVER  (J)  HISTO  CELLS 

419  PRINT  2  LINES  THUS 


422  LET  RMSMC:SQRT.F(RMSMC/REAL.F(NCELLS)) 

423  PRINT  2  LINES  WITH  MAEMC, RMSMC 

424  THUS 

Max  aba  error  and  RMS  error  In  convol  c.d.f.  via  Monte-Carlo  ».»••»»»  »,»»»•»» 

427  'LS'STOP 

428  END  "TEST.CONVOLV 

1  FUNCTION  EFUN  (K,  X) 

2  " 

3  "Test  cum  prob  function  used  In  TEST.CONVOLV.  Function  shown  belci 

4  "  is  a  Erlang  distribution  with  (Integer)  shape  parameter  K  and  stand- 

5  "ardlzed  argument  X. 

6  ■ ' 

7  DEFINE  I,X  AS  INTEGER  VARIABLES 

8  LET  EX:EXP.F(-X) 

9  IF  K=  1 

10  RETURN  WITH  1.0-EX 

1 1  OTHERWISE 

12  LET  FACT: 1.0 


13  LET  XIrl.O 

14  LET  SUMrI.O 

15  FOR  1=1  TO  K-1  DO 

16  LET  FACT=FACT»I 

17  LET  XI=XI»X 

18  ADD  XI/FACT  TO  SUM 

19  LOOP  ’'OVER  I 

20  RETURN  WITH  1.0-EX*SUM 

21  END  "FUNCTION  EFUN 

1  FUNCTION  WFUN  (SHAPE,  ARG) 

2  " 

3  ’’Function  calculates  the  cumulative  probability  function  for  a 

4  "Wei bull  distribution  having  shape  parameter  SHAPE  and  standardized 

5  '* independent  variable  ARG. 

6  " 

7  RETURN  WITH  1.0  -  EXP.F(-ARG»»SHAPE) 

8  END  "FUNCTION  WFUN 


1  FUNCTION  W2WFUN  (ARG) 

2  ’ ' 

3  ' ’Calculates  the  convolution  c.d.f.  with  argument  ARG  of  2  standardized 

4  ’’Weibull  probability  distributions,  each  having  shape  parameter  2. 

5  ” 

6  LET  INTG=EXP.F(0.5«ARG»»2)»(SNORM(ARG).SNORM(-ARG))»ARG« 

7  SQRT.F(PI.C/2,0) 

8  RETURN  WITH  1.0  -  EXP.F(-ARG«»2)*( 1 .0+INTG) 

9  END  "FUNCTION  W2WFUN 

1  ROUTINE  CONVOLV  (NFOLD,  XV,  XV) 

2  ’  ' 

3  "Routine  for  calculating  the  result  of  a  sequence  of  convolutions  on 

4  "2  probability  density  functions  (p.d.f.’s).  The  Ist  density  {XV )  is 

5  "convolved  with  the  2nd  (YV),  and  the  result  is  recursively  convolved 

6  "NFOLD-1  times  with  the  1st  p.d.f.  The  program  returns  the  NFOLD-con- 

7  "convoluted  p.d.f.  (in  complex  form)  in  the  vector  YV.  Method:  The 

8  "program  obtains  the  Fourier  transform  (FT)  of  the  X-series  in  XV(«), 

9  "and  the  Y-series  in  YV(*).  Then,  a  complex  product  is  calculated  and 

10  "placed  in  YV(*).  NFOLD- 1  additional  complex  products  are  taken 

11  "between  YV(*)  and  XV(*).  This  final  product  is  inverted  in  place 

12  "in  YV(»). 

13  '*  NAME  TYPE  ENTRY  VALUE  RETURN  VALUE 

14  "  XV  REAL  ARRAY  COMPLEX  X-DENSITY  FT  OF  X-DENSITY 

15  "  YV  REAL  ARRAY  COMPLEX  Y-DENSITY  CONVO’uUTION  DENSITY 

16  "NOTE;  THE  DIMENSION  Or  ARRAYS  MUST  BE  AN  INTEGER  POWER  OF  2. 

17  " 

18  DEFINE  I,IMAX,K,N,NF0LD,NP2  AS  INTEGER  VARIABLES 

19  DEFINE  XV,  YV  AS  REAL,  1-DIMENSIONAL  ARRAYS 

20  LET  N=DIM.F(XV(»)) 

21  LET  IMAX=N/2 

22  " 

23  "CHECK  VALUE  OF  N. 

24  •  * 

25  LET  NP2=1 

26  'PO'LET  NP2=NP2»2 

27  IF  NP2<N 

28  GO  TO  PO 
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29  OTHERWISE 

30  IF  NP2>N 

31  PRINT  1  LINE  WITH  N  THUS 

IMPROPER  VALUE  Or  INPUT-ARRAY  DIMENSION  {=  »»»••)  IN  ROUTINE  CONVOLV. 

33  STOP 

34  OTHERWISE  "GO  ON 

35  " 

36  "OBTAIN  THE  FOURIER  TRANSFORMS  OF  XV  AND  YV. 

37  " 

38  CALL  F0UR.TRANS(-1,XV(»)) 

39  CALL  F0UR.TRANS(-1,YV(»)) 

40  '  ' 

41  "PLACE  PRODUCT  IN  YV(»). 

42  '  * 

43  FOR  K=1  TO  NFOLD-1  DO 

44  FOR  1=1  TO  IMAX  DO 

43  LET  TEMPR=YV(2»I-1) 

46  LET  TEMPI=YV(2»I) 

47  LET  YV ( 2»I- 1 ) =XV ( 2»I. 1 ) •TEMPR-XV ( 2»I ) •TEMPI 

48  LET  YV ( 2»I ) =  XV ( 2»I- 1 ) •TEMPI+XV ( 2»I ) •TEMPR 

49  LOOP  "OVER  (I)  FOURIER  FREQUENCIES 

50  LOOP  "OVER  (K)  CONVOLUTION  ORDER 

51  " 

52  "GET  INVERSE  TRANSFORM. 

53  " 

54  CALL  FOUR. TRANS (1  ''V:**;; 

55  RETURN 

56  END  "CONVOLV 

1  ROUTINE  FOUR. TRANS  (ISIGN,  DATA) 

2  ' ' 

3  "Routine  to  calculate  the  Fourier  transform  (or  inverse  transform) 

4  "of  a  sampled  data  trace,  which  is  passed  in  the  input  vector  DATA(*). 

5  "The  algorithm  used  is  the  Cooley-Tukey  fast  Fourier  transform  (FFT), 

6  ' ' implemented  by  Norman  Brenner  of  the  MIT  Lincoln  Lab.  The  technique 

7  "requires  that  the  #  of  real  data  points  (N)  be  EXACTLY  2**K,  K  >  0. 

8  "If  ISIGN  5  -1,  the  routine  yields  the  transform.  If  ISIGN  =  1,  the 

9  "inverse  transform  is  produced.  Program  output,  in  either  case,  is 

10  "tne  one-dimensional  array- DATA(*) .  When  giving  the  transform  with 

11  "N/2  complex  frequencies,  requiring  N  elements,  the  real  and  imaginary 

12  "components  are  stored  in  adjacent  storage  positions.  If  a  ISIGN  s  -1 

13  "transform  is  followed  by  a  +1  transform,  the  original  trace  appears 

14  "scaled  by  a  factor  of  N. 

15  "Transform  amplitudes  are  defined  by: 

16  "FT(K)=SUM  OVER  J  :  EXP(.2»PI»IMAG*( J- 1  )»(K.  1  )/N)«DATA( J) , 

17  "1  LE  K  LE  N. 

18  "Input  series  in  DATA  must  be  in  complex  form  with  DIM.F(DATA(*))s2*N. 

19  ’ ' 

20  DEFINE  I,  ISIGN,  NDIM,  AND  N  AS  INTEGER  VARIABLES 

21  DEFINE  IP0,IP1,IP2,IP3,I1,I2A,I2B,I3,AND  I3REV  AS  INTEGER  VARIABLES 

22  DEFINE  DATA  AS  A  REAL,  1-DIMENSIONAL  ARRAY 

2i  LET  NDIM=DTM.F(DATA(»)) 

24  LET  N=MDTM/2 

25  LET  IP0=2 

26  LET  IP3=IP0«N 

27  LET  I3REV=1 

28  FOR  13=1  TO  IP3  BY  IPO  DO  "TO  P50 
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29 
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IF  I3<I3REV 

LET  TEMPR=DATA(I3) 

LET  TEMPI=DATA(I3+1) 

LET  DATA(I3)=DATA(I3REV) 

LET  DATA(I3+1)=DATA(I3REV+1) 

LET  DATA(I3REV)=TEMPR 
LET  DATA(I3REV+1)=TEMPI 
ALWAYS 

LET  IP1=IP3/2 
’P3'IF  I3REV>IP1 

SUBTRACT  IPl  FROM  I3REV 
LET  IP1=IP1/2 
IF  IPl  GE  IPO 
GO  TO  P3 
OTHERWISE 
ALWAYS 

ADD  IP!  TO  I3REV 
LOOP  ’’OVER  13  (P50) 

LET  IP1=IP0 
’P6’IF  IPl  GE  IPS 
RETURN 
OTHERWISE 
LET  IP2=IP1«2 

LET  THETA=2.0»PI. C/REAL. F(ISIGN»IP2/IP0) 
LET  SINTH=SIN.F(THETA/2.0) 

LET  WSTPR=-2.0»SINTH»'»2 
LET  WSTPIzSIN.FCTHETA) 

LET  WRrI.O 
LET  WIrO.O 

FOR  11=1  TO  IP1  BY  IPO  DO 

FOR  13=11  TO  IP3  BY  IP2  DO 
LET  I2A=I3 
LET  I2B=I2A+IP1 

LET  TEMPR=WR»DATA(I2B)-WI»DATA(I2B*1 ) 
LET  TEMPI=WR»DATA ( I2B+ 1 ) ♦WI^DATA ( I2B ) 
LET  DATA(I2B)=DATA(I2A)-TEMPR 
LET  DATA ( I2B> 1 ) =DATA ( I2A^ 1 )-TEMPI 
ADD  TEMPR  TO  DATA(I2A) 

ADD  TEMPI  TO  DATA{I2A^1) 

LOOP  ’’OVER  13 
LET  TEMPRrWR 

LET  WR  =  WR»WSrPR-WI»WSTPUWR 
LET  WI=WI«WSTPR4TEMPR»WSTPI^WI 
LOOP  ’’OVER  II 
LET  IP1=IP2 
GO  TO  P6 

END  ’’FOUR. TRANS 
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